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1.  INTRODUCTION 

1.1  BACKGROUND 

The  development  of  a  missile  guidance  system  generally  requires 
several  phases,  including  inception,  preliminary  design  and  feasibility 
studies,  decisions  concerning  implementation  of  various  system  functions, 
and  compensation  or  design  modification  to  obtain  the  best  possible  system 
performance  under  realistic  constraints.  In  the  later  stages  of  development, 
the  mathematical  system  model  used  as  a  basis  for  generating  system  per¬ 
formance  projections  almost  inevitably  contains  nonlinear  effects  and  random 
disturbances.  Nonlinearity  is  generally  associated  with  nonlinear  relations 
inherent  to  the  laws  of  physics,  unavoidable  hardware  nonlinearities,  and 
essential  design  nonlinearities;  while  random  disturbances  may  include  noise 
(e.g.,  thermal  effects),  sensor  measurement  errors,  random  target  maneu¬ 
vers,  and  random  initial  conditions.  When  random  effects  are  significant, 
some  statistical  measure  of  system  performance  is  required;  often  the  root- 
mean-square  (rras)  miss  distance  achieved  at  the  time  of  target  interception 
is  used  ior  assessing  the  capability  of  a  tactical  missile. 

The  traditional  approach  used  to  obtain  rms  miss  distance  for  guid¬ 
ance  systems  with  significant  nonlinearities  has  been  the  utilization  of  the 
monte  carlo  method.  In  this  technique,  a  large  number  of  computer  simula¬ 
tions  (trials)  are  made  using  the  required  nonlinear  model  with  different 
initial  conditions  and  random  forcing  functions  generated  according  to  the  given 
statistics.  The  resulting  ensemble  of  simulations  provides  the  basis  for  mak¬ 
ing  estimates  of  the  true  rms  miss  distance.  Associated  with  the  monte  carlo 
method  is  the  problem  that  a  large  number  of  trials  is  required  to  provide 
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confidence  in  the  accuracy  of  the  results;  an  ensemble  as  large  as  1000  may 
be  needed  to  obtain  an  accurate  statistical  analysis  for  a  nonlinear  system. 
Thus,  while  the  monte  carlo  method  may  be  useful  for  performing  a  few 
evaluations  of  a  system’s  performance,  it  is  not  a  very  satisfactory  cool  for 
conducting  extensive  sensitivity  and  tradeoff  studies  for  different  values  of 
the  important  guidance  system  parameters,  or  for  conducting  detailed  studies 
of  nonlinear  effects  on  system  performance,  due  to  the  large  expenditure  in 
computer  time  required.  The  purpose  of  this  research  is  to  extend  an  ana¬ 
lytical  technique  for  the  direct  statistical  analysis  of  nonlinear  systems  which 
imposes  a  considerably  smaller  computational  burden.  This  method,  con¬ 
ceived  and  currently  being  developed  by  TASC,  appears  exceptionally  promis¬ 
ing  as  a  means  for  directly  treating  nonlinear  guidance  system  statistical 
behavior;  In  the  sequel  it  is  referred  to  as  the  Covariance  Analysis  DEscrib- 
ing  function  Technique  (CADET™). 

CADET  is  based  on  the  application  of  quasi -linearization  to  permit 
the  application  of  covariance  analysis  to  generate  guidance  system  perform¬ 
ance  statistics.  Thus,  in  treating  nonlinear  systems,  describing  function  theory 
is  used  to  obtain  a  quasi-linear  approximation  for  each  nonlinearity.  The  most 
important  factor  in  the  use  of  quasi-linearization  is  that  one  essential  property 
of  nonlinear  elements  is  retained:  the  amplitude  sensitivity  of  the  input-output 
relation.  Figure  1.1-1  illustrates  the  fundamental  properties  of  the  random 
input  describing  function  for  an  ideal  saturation  or  limiter . 

Previous  work  (Refs.  1  to  3)  has  demonstrated  the  capability  of 
CADET  to  perform  rapid  statistical  analysis  of  complex  nonlinear  guidance, 
systems,  including  the  effects  of  measurement  noise,  parasitic  coupling  due 
to  nonlinear  boresight  error  aberration,  acceleration  command  limiting, 
and  highly  nonlinear  missile  airframe  dynamics.  The  research  being  under¬ 
taken  in  the  present  program  entails  the  generalization  of  the  basic  CADET 
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(•)  NONLINEAR  FUNCTION  (b)  TAYLOR  SERIES  AND  QUASKtNEARIZED  GAINS 

Figure  1.1-1  Illustration  of  Describing  Function  Theory:  the 
Random  Input  Describing  Function  for  a  limiter 

concept  to  render  it  applicable  to  a  very  broad  class  of  missile-target  inter¬ 
cept  problems .  This  will  establish  an  essential  prerequisite  for  the  sys¬ 
tematic  investigation  of  guidance  system  performance  and  design  the 
ability  to  efficiently  perform  detailed  studies  of  the  effects  of  all  significant 
sources  of  system  degradation  and  to  evaluate  guidance  system  modifications 
that  may  correct  or  ameliorate  these  deficiencies .  This  type  of  analysis 
would  be  prohibitively  expensive  in  terms  of  computer  time  without  the  com¬ 
putational  efficacy  inherent  to  CADET. 

While  the  primary  thrust  of  the  development  is  the  extension  and 
refinement  of  an  efficient  tool  for  the  statistical  evaluation  of  the  perform¬ 
ance  of  missile  guidance  systems,  the  overall  scope  of  CADET  is  evidently 
much  more  general.  The  system  model  based  on  a  nonlinear  state  vector 
differential  equation  with  random  inputs  is  of  broad  generality,  being  de¬ 
scriptive  of  many  continuous  systems  with  nondeterministic  disturbances. 
The  specific  nonlinear  effects  studied  herein  are  by  no  means  restricted 
in  occurrence  to  the  missile-target  intercept  problem.  It  is  hoped  that 
the  success  of  the  present  effort  will  encourage  other  applications  of  the 
CADET  concept. 
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1.2  OVERVIEW 

In  order  to  develop  CADET  to  meet  the  need  for  an  efficient  means 
for  generating  performance  projections,  thereby  permitting  the  effective 
diagnosis  and  correction  of  potential  guidance  system  deficiencies,  the 
following  objectives,  have  teen  chosen  for  the  present  study: 

•  Development  of  a  realistic  missile-target  intercept 
model  with  significant  random  effects  and  a  number 
of  quite  diverse  nonlinearities 

•  Development  of  a  monte  carlo  simulation  capability 
to  provide  performance  statistics  for  comparison 
with  CADET  results 

•  Development  of  CADET  methodology  to  permit  the 
efficient  generation  of  performance  statistics  for  the 
given  system  model,  based  on  the  computation  of 
describing  function  approximations  for  each  non¬ 
linearity  and  the  application  of  covariance  analysis 
to  the  quasi-linear  system 

•  Verification  of  the  ability  of  CADET  to  provide  accu¬ 
rate  performance  projections  by  comparing  the  sta¬ 
tistical  analysis  given  by  CADET  with  the  results  of 
corresponding  monte  carlo  analysis 

•  Study  of  the  sensitivity  of  CADET  analysis  to  the 
assumptions  and  approximations  made,  including 
comparisons  of  describing  functions  for  different 
assumed  probability  density  functions  and  compari¬ 
sons  of  histograms  generated  from  monte  carlo 
simulations  with  the  density  functions  assumed  in 
developing  CADET 


The  first  procedure  used  in  developing  confidence  in  the  capability 
of  CADET  to  provide  accurate  performance  statistics  is  a  step-by-step 
study  of  nonlinear  effects.  In  every  phase  of  the  investigation,  the  CADET 
and  monte  carlo  programs  for  statistical  analysis  are  extended  to  include 
identical  system  dynamic  equations,  so  that  performance  projections  are 
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exactly  comparable.  The  comparison  of  CADET  and  monte  carlo  perform¬ 
ance  projections  included  treating  the  following  effects: 

•  The  basic  missile -tar  get  intercept  model,  with  nonlinear 
missile-target  kinematic  relationships,  range  dependent 
error  sources  in  the  seeker  noise  model,  acceleration 
command  limiting,  and  nonlinear  guidance  law 

•  Seeker  mass  imbalance 

•  Nonlinear  friction  in  the  seeker  gimbal 

•  Nonlinear  restoring  torques  acting  on  the  seeker  head 

•  Nonlinear  attenuation  of  boresight  error  due  to  the 
receiver/ signal  processing  nonlinearity 

Following  the  application  of  CADET  to  the  analysis  of  guidance 
system  performance,  several  aspects  of  the  sensitivity  problem  were  con¬ 
sidered: 

•  Sensitivity  of  random  input  describing  function 
calculations  to  changes  in  the  nonlinearity  input 
probability  density  function 

•  Methods  for  calculating  approximate  quasi-linear  gains 

•  Generation  of  approximate  histograms  from  data  pro¬ 
vided  by*  a  large  number  of  monte  carlo  simulations  to 
assess  the  impact  of  deviation  from  the  gaussian  as¬ 
sumption  on  CADET  analysis 

For  reasonable  initial  condition  statistics,  the  CADET-monte  carlo 
comparison  proved  that  C  ADET  results  were  reliable  --  in  fact,  it  appears 
that  the  statistics  given  by  monte  carlo  analysis  are  not  superior  to  the 
CADET  results  until  the  number  of  monte  carlo  trials  is  in  excess  of  several 
hundred.  While  the  sensitivity  of  CADET  to  underlying  assumptions  is 
appreciable  in  some  circumstances,  it  is  gratifying  to  be  able  to  demonstrate 
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that  this  effect  is  quite  small  in  the  present  context  except  in  somewhat  un¬ 
realistic  circumstances  when  nonlinearity  inputs  are  highly  nongaussian. 


1.3  REPORT  OUTLINE 

This  report  is  organized  according  to  the  following  outline:  Chapter  2 
deals  with  the  significant  features  of  the  CADET  and  monte  carlo  methods  for 
generating  performance  statistics  in  nonlinear  systems  with  random  inputs; 
Chapter  3  treats  a  general  discussion  of  the  guidance  system  model;  details 
of  the  initial  verification  procedure  and  case  studies  of  the  subsidiary  non¬ 
linear  effects  are  given  in  Chapter  4;  and  the  analysis  of  the  sensitivity  of 
CADET  calculations  to  incorrect  assumptions  concerning  nonlinearity  input 
probability  density  functions  is  treated  in  Chapter  5.  A  summary  of  the  study 
and  general  conclusions  are  provided  in  Chapter  6.  Appendices  are  included 
to  treat  the  technical  aspects  of  the  system  model  (Appendix  A),  random  in¬ 
put  describing  function  calculation  (Appendix  B),  and  CADET  theory  and 
methodology  (Appendix  C). 
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2.  STATISTICAL  ANALYSIS  VIA  MONTE  CARLO  AND  CADET 


2. 1  THE  MONTE  CARLO  TECHNIQUE 

The  monte  carlo  method  for  the  statistical  analysis  of  the  performance 
of  a  nonlinear  system  with  random  inputs  is  based)  on  direct  simulation  to  deter¬ 
mine  the  system  response  to  "typical”  initial  conditions  and  noise  input  sample 
functions  generated  according  to  their  specified  statistics.  Thus,  the  informa¬ 
tion  required  for  this  analysis  is  the  system  model,  initial  condition  statistics, 
and  random  input  statistics.  The  system  model'  can  be  given  in  the  form  of  a 
state  vector  differential  equation, 

x  =  f  (x)  +  Gw(t)  (2. 1-1) 


where  f(x)  represents  the  linear  and  nonlinear  dynamic  relationships  in  the  sys¬ 
tem,  w  is  an  input  vector,  and  the  matrix  G  specifies  the  input  allocation,  i.  e. , 
each  element*  g.^  of  G  describes  the  effect  of  the  input  element  w^(t)  on  the  state 
variable  derivative  x..  The  state  vector  differential  equation,  Eq,  (2. 1-1),  is 
portrayed  in  block  diagram  notation  in  Fig.  2. 1-1.  The  initial  condition  of  the 
state  vector  is  tpecified  by  assuming  that  the  state  variables  are  jointly  normal 
with  a  given  mean  vector  and  covariance  matrixt, 


E[x(0)]  mQ 
E[x(0)xT(0)]  =  PQ 


(2.1-2) 


*  th  th 

For  any  matrix  G,  the  quantity  gjj  denotes  the  element  in  the  i  row  and  j 

column;  similarly,  for  any  column  vector  w,  w.  denotes  the  element  in  the 
jth  row>  ] 

tE  [  ]  denotes  the  expected  value  of  the  bracketed  variable. 
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Figure  2. 1-1  Nonlinear  System  Model 


The  input  vector  w  is  often  assumed  to  be  composed  of  elements  that  are  inde 
pendent  gaussian  white  noises,  plus  an  additive  mean;  thus 


E[w(t)]  =  b(t) 

E[(w(t)  -  b(t))  (w(t)  -  b(r))T]  =  Q(t)  6(t-r) 


(2. 1-3) 


where  Q(t)  is  the  input  spectral  density  matrix  (which  is  diagonal,  due  to  inde¬ 
pendence)  and  the  impulse  function  5(t  -r)  indicates  that  the  input  vector  random 
components  have  zero  autocorrelation  for  t^r;  i.e.,  the  quantity  u(t)  =  w(t)  -b(t) 
is  white  noise. 

Given  the  above  information,  monte  carlo  analysis  requires  a  large 
number,  say  q,  of  representative  simulations  of  the  system  response,  viz. , 
the  q-fold  repetition  of  the  following  procedure:  First,  an  initial  condition  vec¬ 
tor  is  chosen  according  to  the  statistics  indicated  above;  i.e. ,  a  random  number 
generator  calculates  a  random  vector  x(0)  based  on  Eq.  (2. 1-2).  Then  a  random 
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initial  noise  vector  is  generated,  using  the  statistics  given  in  Eq.  (2. 1-3)*. 
These  vectors  provide  the  data  for  evaluation  of  x(0)  in  Eq.  (2. 1-1)  which  in 
turn  is  used  to  propagate  the  solution  from  t  =  0  to  t  =  h  according  to  any  stan¬ 
dard  technique  for  the  digital  integration  of  a  state  vector  differential  equation. 
Then,  given  x(h),  simulation  continues  by  the  generation  of  a  new  value  of  the 
input  noise  vector  w(h),  evaluation  of  x(h),  numerical  integration  to  obtain  x(2h) 
and  so  on,  to  the  specified  terminal  time  t^ 

Performing  q  independent  simulations  yields  an  ensemble  of  state  tra¬ 
jectories,  each  denoted  x^(t;  x^(0),  w_^(t))  to  stress  the  dependence  of  the 
trajectory  on  the  random  initial  condition  and  noise  input  sample  function: 


x(1)  (t;  x(1)(0),  w(1)(t)) 
x(2)  (t;  x(2\o),  42\t )) 

x<q)  (t;  x^fO),  w(«)(t)) 


0  *  t  S  tf 


(2. 1-4) 


Neglecting  the  bias  component,  we  generate  a  broad -band  gaussian  noise 
uj(t)  with  spectral  density  qj  by  using  a  random  number  generator  to  obtain 
an  independent  sequence  of  gaussian  random  numbers  u.  .,  k  =  0, 1,2, .  .  . 
satisfying  ’ 

E[uijk]  =  0 


E[u1;k]  =qj/h 


Then  we  define  u^(t)  by 


Uj(t)  =  u.  kh  *  t  <  (k  + 1)  h 

where  h  is  a  smalltime  increment.  For  h  small  (1/h  much  larger  than  the 
bandwidth  of  the  system  in  question),  Uj[(t)  is  essentially  a  gaussian  white 
noise  process. 
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Each  satisfies  the  state  vector  differential  equation  (Eq.(2.1-1))  to  within  the 
accuracy  of  the  numerical  integration  method  used.  The  mean  m(t)and  co- 
variance  P(t)  of  the  state  vector  are  estimated  by  averaging  over  the  ensemble 
of  trajectories  using  the  relations 

m(t)  =  i  £  x(i)(t)  a  m(t) 
q  1  =  1 

(2. 1-5) 

P(t)=jrr  £  (x(i)(t)  -  m(t))  (x(1)(t)  -  m(t))T  3  P(t) 
q  1  =  1 

where  m(t)  and  $(t)  denote  the  estimated  values.  The  essence  of  the  monte 
carlo  technique  is  illustrated  in  Fig.  2. 1-2. 


R  *11797 


MONTE  CARLO  SIMULATION  STATE  l:  TRIAL 


Figure  2. 1-2  Schematic  Characterization  of  the 
Monte  Carlo  Technique 
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In  order  to  assess  the  accuracy  of  the  approximate  statistics  given 
in  Eq.  (2. 1-5),  it  is  necessary  to  consider  the  statistical  properties  of  the 
estimates  m(t)  and  P(t).  To  simplify  the  notation,  consider  a  scalar  random 
variable  y  (e.  g. ,  the  value  of  some  system  state  variable  at  some  time  of 
interest),  and  let  m  and  p  represent  the  true  values  of  the  mean  and  variance 

of  y, 


m  =  E  [y] 


p  =  E[(y  -  m)2] 


(2. 1-6) 


By  performing  one  set  of  q  monte  carlo  trials,  we  obtain  a  single  estimate  of 
m  and  p,  which  we  denote  m  and  p.  These  estimates  are  also  random  variables 
that  is,  if  another  set  of  q  monte  carlo  trials  were  performed  independently  of 
the  first  set,  but  with  the  same  statistics  for  the  initial  conditions  and  noise 
inputs,  then  a  different  ensemble  of  simulations  results,  and  different  esti¬ 
mates  for  the  mean  and  variance  would  be  obtained.  If  q  is  sufficiently  large, 

then  we  can  invoke  the  central  limit  theorem  to  justify  the  assumption  that  the 
A  A 

random  variables  m  and  p  are  gaussian*,  and  thus  that  their  distributions  are 
completely  specified  by  the  following  statistics,  asymptotically  true  for  large 
q  and  given  in  Ref.  4: 

r  A, 

E[m]  =  i 
E[p]=I 
'A2  =  E[im-m)2 

=  E[$-p)2] 


n 


I  _  P 

1  "  q 


m4-p 


(2.1-7) 


For  q  <  20,  it  is  necessary  to  assume  that  p  has  the  chi  square  distribution 
if  y  is  a  gaussian  variable  (Ref.  5). 
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I 


P 


where  u.  is  the  fourth  central  moment, 

4 

M4  =  E[(y-m)4]  (2.1-8) 

For  many  common  probability  density  functions  (pdf’s),  the  fourth  central 
moment  is  of  the  form 

P4  =  X  p2  (2. 1-9) 


Table  2. 1-1  gives  a  summary  of  Xfor  the  pdf’s  used  in  this  report.  In  terms 
of  the  parameter  X ,  we  can  express  the  standard  deviations  given  in  Eq. 

(2. 1-7)  in  the  form 


(2.1-10) 


The  above  discussion  of  the  statistics  of  the  gaussian  random  variable 
p  provides  the  basis  for  determining  a  range  in  the  vicinity  of  p  such  that  the 
true  value  of  p  is  guaranteed  to  lie  within  that  range  with  a  specified  probability, 

0.  This  is  done  by  determining  the  number,  n^,  of  standard  deviations,  ap,  such 
that 


Prob  [ 0  ^  |p-p|  *nffcr  ]  =  0 

r 


(2.1-11) 


Since  p  is  a  gaussian  random  variable,  is  the  solution  to 


y/Tv 


~  f  a  exp  (-  i  C2)  dC  =  0 

2#  J_n  2 


(2. 1-12) 
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TABLE  2.  1-1 

SOME  COMMON  PROBABILITY  DENSITY  FACTIONS 


*-1)798 


DESIGNATION 

FUNCTIONAL 
REPRESENTATION  • 

GRAPHICAL 

REPRESENTATION 

X* 

EXPONENTIAL 

-CO  <  X  <  ♦  CO 

.  P 

*) 

6 

m 

NORMAL 

1  /  1  /  x-m 

pv  2  f  cr  ) 

-CO  <*<  ♦  CO 

p 

•/Tna 

_) - »-x 

m 

3 

TRIANGULAR 

s/ba  \  -fba  /' 

tn-JSo  <  x  <  m  *JZa 

P 

*j 

">,X,X<niSS^rn?/^a’ 

2.4 

m-yi(T 

UNIFORM 

_J _ 

/t2ff. 

m-/3<r  <  x  <  m  */5<r 

sfc-j 

(x) 

1.8 

J] 

L. 

m^3tr 

m  rn-Zi^ 

BIPOLAR 

(Oiscralc) 

•j  8  {x-tn-  cr ) 

♦  8  ( x  -m-s  c ) 

j 

*1/2  1 

1.0 

m-<r 

m  m+cr 

•Formulat«d  to  have  mean  m  ond  standard  deviation  a 
t  X  is  referred  to  as  the  "kurlosis  ’*  or  ’’excess"  of  a  pdf.  (Ref.6) 


For  example,  if  the  desired  probability  is  0. 95,  Eq.  (2. 1-12)  yields  =  1. 96. 
Other  values  of  n^  corresponding  to  different  values  of  ip  can  be  obtained  from 
probability  integral  tables  (Ref.  6);  several  representative  values  are  given  in 
Table  2. 1-2. 

TABLE  2. 1-2 

CUMULATIVE  PROBABILITY  WITHIN  n^  STANDARD  DEVIATIONS 

OF  THE  MEAN  FOR  A  GAUSSIAN  RANDOM  VARIABLE 

I - 1 - 1 
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To  reformulate  Eq.  (2. 1-1:2)  into  an  inequality  for  p,  we  make  the  sub¬ 
stitution  for  o  indicated  in  Eq.  (2.  il-10)  into  Eq.  (2. 1-11)  to  obtain  the  equiva- 
P 

lent  statement  that 


Prob 


A 

P 


l+n_  JX-1 


a  V 


P 


(2. 1-13) 


that  is,  the  true  value  of  p  lies  between  the  values  p  and  p  indicated  in  Eq. 

(2. 1-13)  with  probability  0.  Alternatively,  in  terms  of  the  estimated  rms  value 
of  the  variable,  o*,  we  have  the  comparable  result 


Prob  [  a  £  a  *  a  ]  =  0 


for  a  and  cr  given  by 


(2.1-14) 


The  quantities  a  and  o  are  referred  to  as  lower  and  upper  confidence  limits; 
the  value  of  0  expressed  as  a  percent  (1000)  is  the  degree  of  confidence. 
Equation  (2. 1-14)  demonstrates  that  the  standard  deviation  confidence  limits 
can  be  obtained  from  a  simply  by  using  the  multipliers  p  and  p.  The  latter 
are  functions  only  of  the  parameter  X,  the  number  of  monte  carlo  trials  q, 
and  the  number  of  standard  deviations  n  determined  by  the  desired  degree  of 
confidence. 
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The  problem  of  making  a  reasonable  choice  of  X,  which  depends  upon 
the  statistics  of  the  random  variable  y,  must  be  considered  before  the  confi¬ 
dence  limit  multipliers  can  be  calculated.  One  option  is  to  determine  an 
approximate  value  of  X  by  estimating  the  fourth  central  moment  using  the  q 
sample  values  of  the  variable  y,  and  calculating 

A  /A2  a 

X  ^  M4/p  =  X 

The  value  of  X  need  not  be  known  exactly,  since  the  confidence  limits  o  and  a 
are  not  very  sensitive  to  errors  in  this  parameter.  In  the  absence  of  reliable 
information  about  the  higher  central  moments,  it  is  frequently  assumed  that  y  is 
gaussian;  i.e.,  that  x  =  3.  However,  if  there  is  reason  to  believe  that  the  pdf 
for  y  has  abnormally  heavily  weighted  tails  —  as  in  the  case  of  the  exponential 
distribution  in  Table  2. 1-1  --  then  a  larger  value  of  X  may  be  required. 

Typical  values  of  p  and  p  ior  X  =  3  are  indicated  as  functions  of  the 
number  of  monte  carlo  trials  in  Fig.  2. 1-3,  for  two  values  of  confidence.  As 
an  example  of  the  significance  of  the  confidence  interval,  if  we  desire  to  have 
99%  certainty  that  a  is  within  10%  of  the  estimated  value,  o;  i.e. , 

Prob  [0.90  o  *  a  *  1.1  a]  =  0.99  (2.1-15) 

then  Fig.  2.1-3  demonstrates  that  it  is  necessary  to  perform  440  trials; 

256  trials  suffice  for  95%  confidence.* 

Figure  2. 1-4  shows  the  deterioration  that  occurs  in  the  accuracy  of 
the  monte  carlo  estimated  standard  deviation  for  a  given  level  of  confidence  if 
X  is  greater  than  three  due  to  y  being  nongaussian.  We  discuss  an  instance 
where  X  =  19.  5  in  Section  5. 3;  in  this  case,  even  for  more  than  200  trials,  the 
upper  95%  confidence  limit  is  nearly  50%  greater  than  the  estimated  value  of  a. 


$ 

Note  that  the  bounds,  pand  p,  are  not  symmetric  with  respect  to  one; 
thus  the  point  at  whictTp  crosses  1. 1  determines  the  value  of  q  for  which 
Eq.  (2.1-15)  is  satisfied. 
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Figure  2. 1-3  Typical  Confidence  Interval  Multipliers 
for  the  Estimated  Standard  Deviation  of  a 
Gaussian  Random  Variable  (X  =  3) 


The  confidence  interval  calculation  for  the  estimated  mean  is  quite 
direct,  since  cr/\  (Eq.  (2. 1-9))  is  not  a  function  of  the  mean.  The  same  value 
of  n^  is  obtained  for  the  desired  degree  of  confidence  (e.g. ,  n^  =  1. 96  for  95% 
confidence),  and  value  of  p  given  in  Eq.  (2. 1-13)  is  used  in  deriving  the  result 
that 


Prob  [  m  ^  m  £  m  ]  =  0 
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The  confidence  limit  concept  developed  above  provides  a  statistical 
measure  of  the  accuracy  of  the  estimated  mean  and  standard  deviation  of  a 
random  variable  obtained  by  using  the  monte  carlo  method.  It  is  only  possible 
to  assess  the  accuracy  of  such  estimates  in  a  probabilistic  sense;  e.  g. ,  for 
256  trials,  we  can  assert  that  an  estimated  standard  deviation  (rms  value)  is 
within  10%  of  the  true  value,  with  probability  0.95  (with  95%  confidence). 

2. 2  THE  COVARIANCE  ANALYSIS  DESCRIBING  FUNCTION  TECHNIQUE 
(CADET) 

Covariance  analysis,  where  it  is  applicable,  provides  a  direct,  exact 
technique  for  the  statistical  evaluation  of  the  performance  of  dynamic  systems 
with  random  inputs,  permitting  the  propagation  of  the  mean  component  and  co- 
variance  matrix  of  the  system  state  vector  as  functions  of  time.  This  technique 
does  not  require  the  generation  of  a  large  ensemble  of  representative  state  tra¬ 
jectories  and  the  computation  of  the  ensemble  statistics.  The  latter  approach 
—the  monte  carlo  method  —  is  both  time  consuming  (in  terms  of  computer 
time)  and  approximate;  we  have  seen  in  Section  2. 1  that  several  hundred  sample 
state  trajectories  may  be  required  in  order  to  obtain  reasonably  accurate  sta¬ 
tistics  for  systems  in  which  the  state  variables  are  nearly  gaussian;  for  non¬ 
linear  systems  in  which  variables  may  be  highly  nongaussian,  one  thousand  or 
more  simulations  could  be  required  to  achieve  an  acceptable  level  of  confidence 
in  the  accuracy  of  the  results.  The  direct  approach  is  thus  distinctly  preferable 
to  the  monte  carlo  technique. 

The  fundamental  bases  of  covariance  analysis  are  the  differential  equa¬ 
tions  governing  the  evolution  of  the  mean  vector  and  covariance  matrix  with 
time.  Corresponding  to  the  state  vector  differential  equation  given  by  Eq.  (2. 1-1) 
with  the  input  vector  specified  by  Eq.  (2. 1-3),  we  have  the  differential  equations 
(Ref.  7) 
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0 


1 5r 


0 


n 


U 


m  =  E[f(x)]  +  Gb(t) 

P  =  E [f  rT]  +  E[r  fT]  +  GQ(t)G 


where 


E[x(t)]  =  m(t) 

r(t)  =  x(t)  -  m(t)  (2. 2-2) 

P(t)  =  E[r(t)  rT(t)] 


The  first  and  second  moments  of  the  system  response  are  completely  deter¬ 
mined  by  the  integration  of  the  indicated  vector  and  matrix  differential  equations 
(Eq.  (2.2-1))  when  the  initial  conditions,  m(0)  and  P(0),  are  specified. 

The  form  of  the  differential  equations  for  the  statistics  of  x(t)  is  par¬ 
ticularly  convenient  in  the  case  of  linear  systems;  given 

f(x)  =  Fx 


Eq.  (2.2-1)  reduces  to 


m  =  Fm  +  Gb 
P  =  FP  +  PFT  +  GQGT 


(2.2-3) 


These  linear  differential  equations  can  be  solved  using  standard  numerical  inte 
gration  techniques. 
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The  nonlinear  covariance  equations  can  be  represented  in  the  same  for¬ 
mat  as  Eq.  (2. 2-3)  by  definition  of  two  auxiliary  matrices,  and  Nr,  given  by 


Nm21  =  E[f(x)] 

NrP  =  E[f(x)  rT] 

so  Eq.  (2. 2-1)  may  be  written 

m 

P 

=  2  sym  [NrP  ]  +  Q 


N  m  +  Gb 

ni —  — 


=  N  P  +  PNT  +  GQGT 


(2. 2-4) 


(2.2-5) 


where  sym  [  ]  indicates  the  symmetric  part  of  the  indicated  matrix,  and  Q  is 
the  spectral  density  matrix  of  G(w  -  b).  The  dependence  of  Nm  and  on  the 
statistics  of  the  state  vector,  which  is  implicit  in  the  expectation  operations  in 
Eq.  (2. 2-4),  is  due  to  the  existence  of  nonlinearities  in  the  system.  Without  non¬ 
linear  effects,  the  propagation  of  the  mean  and  covariance  is  "uncoupled, "  as 
shown  in  Eq.  (2. 2-3). 

The  matrices  N  and  N  given  in  Eq.  (2. 2-4)  must  be  determined  be- 
m  r 

fore  we  can  proceed  to  solve  Eq.  (2. 2-5).  A  direct  approximate  method  of 
accomplishing  this  is  to  use  describing  function  theory  to  linearize  the  system 
nonlinearities;  the  resulting  generalization  of  linear  covariance  analysis  is  then 
called  CADET  —  Covariance  Analysis  Describing  function  Technique.  This  pro¬ 
cedure,  presented  in  detail  in  Appendix  C,  entails  assuming  the  form  of  the  non¬ 
linearity  input  and  calculating  describing  functions  which  provide  a  quasi -linear 
approximate  input-output  relation  for  each  nonlinearity.  For  the  problem  at  hand, 
the  system  input  vector  w(t)  is  taken  to  be  a  gaussian  random  process  plus  a 
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bias,  and  it  is  useful  to  assume  that  the  state  variables  also  have  mean  com¬ 
ponents  and  random  parts  that  are  jointly  normal.  While  this  assumption  is 
strictly  true  only  for  linear  systems,  it  is  often  approximately  valid  in  non¬ 
linear  systems.  Although  the  output  of  a  nonlinearity  with  a  gaussian  input  is 
generally  nongaussian,  it  is  known  from  the  central  limit  theorem  that  random 
processes  tend  to  be  made  gaussian  when  passed  through  low-pass  linear  sys¬ 
tems  ("filtered").  Hence  if  there  are  a  few  stages  of  linear  dynamics  between 
nonlinearities,  the  input  to  each  nonlinearity  should  be  nearly  gaussian.  The 
essential  requirement  of  CADET  is  that  the  state  variables  must  be  nearly 
jointly  normal.  From  a  practical  viewpoint,  the  gaussian  hypothesis  serves  to 
simplify  the  CADET  methodology  by  permitting  each  nonlinearity  to  be  treated 
in  isolation,  with  N  and  N  formed  from  the  individual  random  input  describ- 
ing  functions  (ridf's)  for  each  nonlinearity,  as  discussed  in  Appendix  C.  Since 
ridf’s  have  been  extensively  catalogued  in  Ref.  8,  the  implementation  of  CADET 
is  a  straightforward  procedure  for  a  broad  class  of  nonlinear  systems.  We  also 
note  that  under  the  gaussian  assumption,  the  random  input  describing  functions 
can  be  calculated  directly  from  the  mean  vector,  m,  and  the  covariance  matrix, 
P,  of  the  system  state  vector.  Thus,  we  write  N  and  N  in  the  form 


=  N  (m,  P) 
m  — ’  ' 


Nr  =  Nr(m,  P) 


Relations  of  this  form  permit  the  direct  evaluation  of  the  ridf’s  at  each  integra¬ 
tion  step  in  the  propagation  of  m  and  P  by  numerical  techniques  according  to 
Eq.  (2.2-5),  as  illustrated  in  Fig.  2.2-1. 

To  demonstrate  the  ease  with  which  CADET  can  be  mechanized  under 
the  gaussian  assumption,  we  consider  a  single  nonlinearity  and  illustrate  the 
steps  involved  in  performing  statistical  analysis  via  CADET.  Assume  that  only 
the  nonlinearity  f.  (x. )  occurs  in  the  differential  equation  for  x., 
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R- 11 805 


Figure  2.  2-1  Nonlinear  Covariance  Analysis  --  CADET 
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• 

=  Fx  + 

• 

+  Gw(t) 

(2.2-6) 

• 

• 

• 

X 

n 

• 

0 

where  we  have  explicitly  separated  the  linear  portion  of  the  dynamics  from  the 

nonlinearity.  The  random  input  describing  function  approximation  to  f .(x^)  is 

shown  schematically  in  Fig.  2.  2-2,  The  input  is  first  separated  into  its  mean 

and  random  components,  mfc  and  r^,  which  are  then  multiplied  by  the  scalar 

quasi-linear  gains,  nfk  and  nmk;  the  latter  are  selected  to  minimize  the  mean 

square  error  between  the  linear  approximation  nrkrk  +  nmkmk  and  i(xk).  The 

gain  n  .  is  added  to  the  element  f,,  of  the  matrix  F  in  determining  N  ;  similarly 
rK  jk  r 


*  T 

The  indicated  symmetry  operator  yields  2  sym[N  P]  =  NrP+PN  where  the 

symmetry  of  P  is  taken  into  account.  1 
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nrkrk  +  "mkmk£*j(xkl 


Figure  2.2-2  Random  Input  Describing  Function 
Approximation  of  fjtefc) 

n  ,  is  added  to  the  same  element  of  F  to  give  N  .  Thus,  tne  jth  rows  of  N 
mk  m  r 

and  Nm  are  given  by 


[Nr]j  ~  [fji>  fj2»  *  *  *  ’  (fjk  +  nrk)j  *  *  • '  fj J 

[NiJj  =  tfjl’  V  *  *  *  ’  (fjk  +  nmk^  *  *  * ’ fj  J 


A  comparison  of  quasi -linearization  with  the  classical  Taylor  series  or 
small-signal  linearization  technique  provides  a  great  deal  of  insight  into  the 
success  of  the  ridf  in  capturing  the  essence  of  nonlinear  effects.  Small-signal 
linearization  for  a  scalar  nonlinear  element  f(x)  is  based  on  the  identification  of 
a  nominal  operating  point  xQ  and  the  evaluation  of  the  slope  of  the  nonlinearity 
at  that  value;  then  the  approximation  is  made  that 

f(x)  =  f(xQ)  +  f'(xo)(x-xo)  (2.2-7) 

which  represents  the  first  two  terms  of  a  Taylor  series  expansion  about  the 

3 

given  operating  point,  as  illustrated  in  Fig.  2.2-3  for  the  example,  y=x  .  While 
this  is  a  useful  approach  if  excursions  from  the  nominal  are  small,  the  validity 
of  the  Taylor  series  approximation  Is  questionable  when  x  is  a  random  variable 
which  can  exhibit  large  variations  about  its  mean  value. 
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Figure  2.2-3 


Taylor  Series  Linearization  of  y  =  x 


3 


about  xQ  =  l 


By  contrast,  the  quasi -linear  representation  of  a  nonlinearity  is  sensi¬ 
tive  to  the  input  amplitude  in  some  sense;  in  the  case  of  random  inputs,  the 

2 

statistics  m  =  E[x]  and  p  =  E[(x-m  )  ]  provide  the  measure  of  input  ampli- 
tude.  For  the  example  y =x**,  where  x  is  a  gaussian  random  process,  we  calcu¬ 
late  the  describing  functions  according  to  Eqs.  (C.  2-4)  and  (C.  2-5), 

n  =  3(p  +  m  2) 
r  *x  x 


n 


m 


=  3p  +  m 
x  x 


and  the  nonlinearity  is  approximated  by 

x3  =  (3px  +  mx2)  mx  +  3(px  +  m^)  (x - mx)  (2. 2-8) 

Comparing  Eqs.  (2. 2-7)  and  (2. 2.  8),  we  see  that  the  describing  function  gains 
depend  on  both  the  mean  and  variance  of  x,  as  indicated  in  Fig.  2. 2-4,  while  the 
coefficients  in  the  Taylor  series  approximation  do  not. 

Finally,  some  comments  about  the  generality  of  CADET  are  in  order. 
Many  approximate  techniques  for  treating  nonlinear  systems  are  applicable  only 
to  low  order  systems  (with  two  or  three  state  variables  at  most)  with  one 
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(a)  MEAN  COMPONENT  GAIN  (b)  RANDOM  COMPONENT  GAIN 

3 

Figure  2. 2-4  Quasi-Linearization  of  y  =  x  for  Unity  Input  Mean 

nonlinearity.  By  contrast,  describing  function  techniques  can  be  used  for  sys¬ 
tems  with  any  number  of  nonlinearities,  leading  to  a  quasi -linear  system  model 
that  can  then  be  handled  using  the  CADET  methodology.  There  is  no  limit  on 
the  number  of  state  variables  in  the  system  model;  in  fact,  CADET  tends  to  be 
more  accurate  as  the  ratio  of  the  number  of  states  to  nonlinear  elements  in¬ 
creases,  thus  improving  the  validity  of  the  gaussian  assumption. 

2. 3  COMPARISONS  AND  PHILOSOPHY  OF  APPLICATION 

In  comparing  CADET  and  monte  carlo  methods  for  use  in  obtaining 
performance  projections  for  nonlinear  systems  with  random  inputs,  there  are 
several  significant  similarities.  Both  techniques  are  applicable  to  nonlinear 
system  models  with  an  arbitrary  number  of  states  and  nonlinearities,  and  we 
often  rely  on  the  gaussian  assumption  in  assessing  the  accuracy  of  the  perfor¬ 
mance  statistics  obtained.  In  both  cases,  any  departure  from  normality  can  be 
compensated  for  to  a  certain  extent;  in  CADET,  nongaussian  pdf's  can  be  used 
in  calculating  describing  functions,  while  in  monte  carlo  simulation  the  fact 
that  the  confidence  limits  increase  for  nongaussian  random  variables  (Fig. 

2. 1-4)  can  be  counteracted  by  increasing  the  number  of  trials  performed.  The 
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principal  trade-off  between  the  two  methods  is  in  efficiency  versus  generality 
(versatility).  The  monte  carlo  simulation  ensemble  of  q  representative  state 
trajectories  ^Eq.  (2. 1-2))  can  be  used  not  only  as  a  data  base  for  calculating 
estimated  performance  statistics  m(t)  and  P(t)  at  instants  of  time  of  interest, 
but  also  for  estimating  higher  order  moments,  and  for  generating  histograms 
which  are  approximate  pdf's  for  the  variables  under  consideration  (of  course, 
at  an  additional  cost  in  terms  of  data  processing).  CADET,  on  the  other  hand, 
provides  approximate  values  for  m(t)  and  P(t)  in  a  single  numerical  integration 
of  the  quasi-linear  covariance  equations  ^Eq.  (2. 2-5)),  usually  in  a  small  frac¬ 
tion  of  the  computer  processing  time  required  for  accurate  monte  carlo  analysis. 
Depending  on  the  number  of  state  variables  and  nonlinearities,  and  on  the  de¬ 
sired  accuracy  (which  determines  the  number  of  monte  carlo  trials  required), 
it  may  be  possible  to  perform  from  ten  to  thirty  distinct  CADET  studies  at  the 
same  computational  expense  required  for  one  monte  carlo  analysis. 

One  of  the  primary  purposes  of  the  statistical  analysis  ot  tactical 
missile  system  performance  is  the  evaluation  of  guidance  effectiveness  with 
variations  in  random  input  levels,  Initial  conditions  statistics,  system  parame¬ 
ter  values  and  secondary  nonlinear  effects  such  as  seeker  mass  imbalance, 
acceleration  command  limiting,  etc.  Due  to  the  multiplicity  of  these  factors,  it 
is  evident  that  the  analysis  will  generally  be  done  repeatedly.  As  a  consequence, 
efficiency  is  an  important  consideration;  this  point  is  a  strong  argument  in  fa¬ 
vor  of  CADET.  On  the  other  hand,  the  versatility  of  monte  carlo  simulation 
provides  a  self -check  capability;  i.  e. ,  the  q  representative  state  trajectories 
can  be  used  to  estimate  higher  moments  and  the  pdf’s  of  the  state  variables, 
which  in  turn  leads  to  an  assessment  of  the  accuracy  of  the  monte  carlo  analysis. 
This  is  a  feature  lacking  in  CADET  which  makes  it  advisable  to  rely  upon  monte 
carlo  simulation  in  a  monitoring  capacity,  since  it  is  always  possible  to  obtain 
accurate  performance  projections  by  increasing  q  sufficiently. 
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The  effective  use  of  CADET  and  monte  carlo  analysis  in  concert  can 
be  demonstrated  in  a  hypothetical  trade-off  study  where  two  parameters,  say 
aj  and  0^  to  be  varied  over  certain  ranges  to  obtain  optimal  performance 
in  some  sense  (perhaps  to  obtain  minimum  rms  terminal  miss  distance).  As 
shown  in  Fig.  2. 3-1,  a  few  points  in  the  parameter  plane  are  chosen  for  careful 
CADET -monte  carlo  comparison  (verification  of  CADET);  then  extensive  per¬ 
formance  curves  are  generated  using  CADET,  from  which  the  optimal  values  of 
Cj  and  cCg  are  chosen.  If  desired,  the  vicinity  of  the  point  of  optimally  can  be 
studied  using  a  few  selected  values  of  and  and  performing  the  required 
monte  carlo  simulations.  Similar  approaches  can  be  used  in  studying  sensi¬ 
tivity  to  nonlinear  and  random  effects. 


Figure  2. 3-1  Illustration  of  CADET  and  Monte  Carlo 
Analysis  in  a  Parameter  Trade-Off  Study 
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The  overall  philosophy  of  CADET  usage,  based  on  the  strong  points  of 
both  CADET  and  monte  carlo  simulation,  is  illustrated  in  Fig.  2.3-2.  The 
initial  verification  procedure  is  generally  undertaken  for  the  "nominal  system," 
i.  e. ,  for  the  system  with  nominal  parameters  values,  and  is  of  necessity  more 
meticulous.  Thus  several  hundred  monte  carlo  trials  may  be  performed,  and 
if  there  is  reason  to  believe  that  the  system  is  highly  nonlinear  —  so  that  the 
system  variables  may  be  quite  nongaussian  —  it  may  be  necessary  to  investi¬ 
gate  higher  order  moments  or  histograms  to  decide  whether  more  trials  are 
needed  in  order  to  obtain  a  reliable  statistical  analysis.  Once  this  phase  has 
been  completed  satisfactorily,  the  CADET  parameter  sensitivity  studies  can 
then  be  performed.  Observe  that  the  preliminary  careful  but  time-consuming 
monte  carlo  study  (perhaps  with  the  indicated  self -check  procedure  and  increase 
in  the  number  of  trials)  is  always  required  if  accurate  performance  statistics 
are  to  be  obtained  from  monte  carlo  simulation  with  high  confidence.  The  sub¬ 
sequent  use  of  the  monte  carlo  statistics  to  verify  CADET  requires  minimal 
computer  time  (if  any),  and  that  comparison  paves  the  way  for  the  ensuing 
efficient  study  of  various  effects  via  CADET.  In  the  latter  studies,  it  is  rarely 
necessary  to  perform  as  detailed  a  monte  carlo  analysis  as  is  required  in  the 
initial  verification  procedure,  so  as  few  as  20  to  50  monte  carlo  trials  may 
suffice  to  demonstrate  that  CADET  has  accurately  captured  the  effect  in  ques¬ 
tion.  This  approach  mirrors  the  development  of  CADET  that  has  been  carried 
out  in  the  present  study. 
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Figure  2.3-2  Philosophy  of  CADET  Application 
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3.  '  MISSILE -TARGET  INTERCEPT  MODEL  DEVELOPMENT 


3. 1  INTRODUCTION 

Previous  development  of  CADET  as  an  efficient  tool  for  undertaking 
the  statistical  analysis  of  tactical  missile  guidance  system  performance  lias 
focused  on  the  planar  missile -target  intercept  problem.  At  each  stage,  the 
model  has  been  made  more  general  by  the  inclusion  of  more  of  the  system 
dynamics  (addition  of  more  system  states)  and  the  consideration  of  more  non¬ 
linear  effects  and  random  disturbances.  The  initial  work  (Ref.  1)  treated  only 
one  important  nonlinearity  —  acceleration  command  limiting  —  and  one  noise 
input  —  a  random  target  maneuver.  A  second,  more  extensive  application  of 
CADET  (Ref.  2)  included  three  nonlinear  effects  and  four  random  inputs,  two 
of  which  had  a  deterministically  time -varying  rms  value.  A  third  investiga¬ 
tion  (Ref.  3),  which  complements  and  in  some  respects  parallels  the  present 
effort,  involves  a  very  detailed  application  of  CADET  to  a  specific  tactical 
missile  —  the  SAM-D.  A  brief  synopsis  of  the  system  models  considered  in 
the  various  research  programs  is  provided  in  Table  3. 1-1. 

The  purpose  of  the  present  study  is  to  significantly  extend  the  verified 
capability  of  CADET  to  provide  rapid,  accurate  assessments  of  tactical  missile 
performance.  We  have  generalized  the  previous  missile-target  intercept  models 
by  adding  a  variety  of  nonlinear  effects  and  random  inputs  which  have  heretofore 
not  been  treated  and  which  can  have  a  significant  impact  on  the  effectiveness  of 
the  tactical  missile  —  as  quantified  by  the  rms  miss  distance  between  the  mis¬ 
sile  and  target  at  terminal  time.  The  specific  equations  that  represent  the 
dynamics  of  the  missile  and  target  are  developed  in  Appendix  A;  in  this  chapter, 
we  provide  a  brief  summary  of  the  effects  studied.  The  verification  that  CADET 
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can  indeed  provide  accurate  performance  projections  with  these  effects  taken 
into  account  is  presented  in  Chapters  4  and  5. 

3.2  MISSILE -TARGET  KINEMATICS 

hi  Fig.  3. 2-1  we  portray  the  graphical  definition  of  the  coordinate 
frame  and  variables  involved  in  describing  the  motion  ot  the  missile  and  target. 
In  deriving  the  equations  of  motion,  it  is  assumed  that  the  missile  and  target 
velocity  vector  magnitudes  are  constant,  or,  equivalently,  that  the  missile  and 
target  acceleration  vectors  are  normal  to  the  velocity  vectors.  These  condi¬ 
tions,  neglecting  the  effects  of  drag  and  assuming  that  the  angle  of  attack  is 
small,  are  representative  of  many  missile-target  engagement  situations  during 
the  critical  last  few  seconds.*  Thus,  the  lateral  acceleration  of  either  vehicle 
simply  produces  a  rotation  of  the  corresponding  velocity  vector,  given  by 


y -AXIS  R-11592 


Figure  3. 2-1  Missile -Target  Planar  Intercept  Geometry 


7 

The  application  of  CADET  to  cases  where  angle  of  attack  and  drag  variations 
are  important  is  currently  being  considered  in  a  continuation  of  the  study 
described  in  Ref.  3. 
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(3.2-1) 


The  equations  describing  the  motion  of  the  target  with  respect  to  the  missile 
center  of  gravity  (governing  the  cross-  and  down-range  missile -target  separa¬ 
tion,  x  and  y  respectively)  are  found  by  projecting  the  velocity  vectors  onto 
the  axes  shown  in  Fig.  3. 2-1;  in  terms  of  the  velocity  magnitudes  v  and  v, , 

in  t 


X=-Vm  COS  (e£)  -  Vt  COS  (ea) 
y  =  -  vm  sin  (6^)  +  vt  sin  (e&) 


(3.2-2) 


Equation  (3. 2-2)  represents  the  essential  nonlinearities  inherent  to  the  missile - 
target  kinematic  relationship. 


3.  3  TARGET  LATERAL  ACCELERATION  MODEL 


An  important  source  of  missile  guidance  system  error  is  target  maneu¬ 
verability.  We  assume  that  the  target,  lateral  acceleration  magnitude,  a^,  is  a 
first-order  Markov  process,  modelled  as  a  zero-mean  gaussian  white  noise* 
passed  through  a  single  stage  of  low -pass  filtering,  as  depicted  in  Fig.  3. 3-1.  t 


The  five  white  noise  inputs  to  the  system  are  simply  denoted  w.,  j  =  1, 2, . . . ,  5, 
for  convenient  reference.  ^ 

f Although  the  transfer  function  representation  is  formally  restricted  to  use  in 
completely  linear  systems,  we  use  it  as  a  compact  notation  for  depicting  linear 
subsystems  in  our  model.  While  general  practice  is  to  replace  the  variables 
W5  and  a^  with  the  formal  Laplace  transforms,  say  W^s)  and  At(s),  we  con¬ 
tinue  to  indicate  linear  subsystems  inputs  and  outputs  by  their  time  domain 
representation  to  avoid  a  cumbersome  dual  notation  for  all  variables. 
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(a)  DIFFERENTIAL  EQUATION 
REPRESENTATION 
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(b)  TRANSFER  FUNCTION 
FORMULATION 

Figure  3, 3-1  Band-Limited  Gaussian  Noise 
Model  for  Target  Lateral 
Acceleration 


By  adjusting  the  values  of  target  maneuver  bandwidth,  and  rms 
level,  a  wide  range  of  target  maneuver  characteristics  can  be  represented.  A 
constant  rms  level  of  lateral  acceleration  over  the  entire  terminal  homing 
phase  is  assumed, 

E[at3(t)]  =  <rt2  (3.3-1) 

Thus  we  choose  the  initial  condition  and  gaussian  white  noise  statistics  to 
satisfy 

E[at2(0)]  =  at2 

(3.3-2) 

E[w5(t)wg(T)]  =  2(ot<Tt2  fi(t-T)  =  qg  6  (t  -r) 
where  q,.  is  the  spectral  density  of  this  white  noise  input. 


3-5 


THE  ANALYTIC  SCIENCES  CORPORATION 


3. 4  AUTOPILOT  -  AIRFRAME  MODEL 

We  use  a  linear  time -invariant  autopilot -airframe  model,  under  the 
assumptions  that:  * 

•  Missile  velocity  is  constant  (drag  effects  are 
negligible  over  the  period  of  time  considered). 

•  Altitude  remains  nearly  constant. 

•  The  center  of  pressure,  mass  and  inertia  of 
the  missile  are  constant. 

•  Lift  force  is  linearly  related  to  changes  in  angle  of 
attack  about  some  trim  condition  and  to  control 
fin  deflection. 

•  Control  fin  actuator  dynamics  are  linear. 

•  Fin  effectiveness  is  independent  of  angle  of 
attack. 

The  airframe  equations  of  motion  are  then  based  on  a  set  of  missile  aerody¬ 
namic  coefficients  identified  with  the  vehicle  short -period  dynamics  with 
values  typical  of  a  tactical  missile  in  the  terminal  homing  phase.  We  choose 
the  autopilot  feedback  compensation  gains  to  achieve  a  suitable  compensated 
airframe  response  (Appendix  A,  Section  A.  4).  The  outputs  of  the  airframe 

model  are  missile  lateral  acceleration,  am,  and  missile  body  angular  rate, 

• 

0m;  in  transfer  function  form,  they  are  related  to  the  acceleration  command, 

a  ,  provided  by  the  guidance  law,  by 
c 

720  -  0.  865  s  -  1.  87  s2 
720  +  275  s  +  18.3  s2  +  s3  c 

0.24  +  0.642  s 

720  +  275  s  +  18.3  s2  +  s3  c 

Application  of  CADET  to  a  complex,  nonlinear  airframe  model  is  being  con¬ 
sidered  in  a  continuation  of  the  work  described  in  Ref.  3. 


a 


m 


m 
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which  corresponds  to  a  dominant  pole  at  s  =  -  3. 16  and  secondary  poles  at 
s  =  -  7.  56  ±  13.  Oj. 


3.  5  GUIDANCE  LAW  MODEL 


The  present  study  utilizes  a  nonlinear  model  of  the  classical  propor¬ 
tional  guidance  law,  described  in  Section  A.  5.  In  principle,  it  is  desired  to 
develop  a  missile  lateral  acceleration  that  has  a  component  normal  to  the 
line-of-sight  (LOS)  which  is  proportional  to  the  product  of  closing  velocity  and 
LOS  angular  rate.  In  practice,  it  is  not  possible  to  obtain  either  the  closing 
velocity  or  the  LOS  angular  rate  exactly  for  generating  an  acceleration  com¬ 
mand.  A  noisy  measurement  of  LOS  angular  rate,  denoted  by  v  in  Fig.  3. 5-1, 
is  provided  by  the  missile  seeker,  which  is  treated  in  Section  3. 6.  In  the  gui¬ 
dance  law  model,  a  single  stage  low-pass  filter  is  included  to  reduce  the  effect 

A 

of  measurement  noise  in  rj;  the  filter  output  is  denoted  by  0.  Measurements 
of  closing  velocity  (range  rate)  will  also  contain  errors,  which  are  repre¬ 
sented  by  the  uncertainty  variable  ev-  The  latter  is  a  first-order  markov 
process,  *  modelled  in  the  same  form  shown  in  Fig.  3.  3-1.  The  resulting 
model  for  the  acceleration  command  is  given  by 


cos(0  +0). 
+  vt  cos(0*-0). 


(3.5-1) 


where  the  constant  n;  is  designated  the  navigation  ratio.  The  guidance  law  is 
completed  by  noting  that  this  unconstrained  acceleration  command  must  gen¬ 
erally  be  limited  in  magnitude  in  order  to  prevent  exceeding  the  structural 
capacity  of  the  airframe  and  to  avoid  the  possibility  that  the  missile  might  go 


A  constant  or  bias  error,  eVb,  is  obtained  by  suitably  choosing  the  initial 
condition,  ey(0)  =  eVfo,  and  mean  value  of  the  input,  b 4  =  E[w4j  =  <04  e^. 
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into  a  stall.  Thus  the  actual  acceleration  command,  ac,  is  represented  as 
the  output  of  a  limiter  whose  input  is  a 7  . 


CLOSING  VELOCITY 
ERROR  MODEL 

f - * - , 


1-115990. 


NOISE  FILTER  CLOSING  VELOCITY  WITH  SECANT  ACCELERATION 

COMPENSATION  AND  UNCERTAINTY  COMMANDS 


Figure  3.  5-1  Guidance  Law  Model 
3.  6  THE  SEEKER  SUBSYSTEM 

The  basic  function  of  the  seeker  is  to  track  the  target  and  provide  the 
guidance  package  with  an  estimate  of  the  LOS  angular  rate.  There  are  a  sig¬ 
nificant  number  of  important  error  sources  which  can  lead  to  a  marked  dete¬ 
rioration  in  the  seeker  performance.  We  can  categorize  them  as  noise  sources, 
boresight  error  distortion  effects,  and  seeker  head  disturbance  torques;  there 
are  several  specific  mechanisms  that  give  rise  to  errors  of  each  type. 

•  Noise  sources: 

Target  amplitude  and  angular  scintillation  noise 
Receiver  noise 
Seeker  servo  noise 
External  jamming 
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•  Boresight  error  distortion  sources: 

Aberration 

Receiver/ signal  processing  characteristics 

•  Disturbance  torque  sources: 

Seeker  mass  imbalance 

Seeker  gimbal  friction 

Seeker  head  restoring  torques 

The  seeker  noise  model  includes  three  fundamental  types  of  additive 
random  inputs:  those  that  have  an  effective  rms  level  that  varies  directly  with 
the  missile -target  range,  or  range  proportional  noise,  including  such  effects 
as  receiver  noise  and  distant  external  stand-off  jamming;  those  exhibiting  a 
constant  effective  rms  level,  or  range  independent  noise,  such  as  target  am¬ 
plitude  scintillation  and  seeker  servo  noise;  and  those  with  an  effective  rms 
level  that  varies  inversely  with  the  range,  or  inverse  range  proportional  noise, 
such  as  target  angular  scintillation.  These  sources  are  modelled  as  wide-band 
independent  gaussian  noises  (filtered  gaussian  white  noise)  with  appropriate 
nonlinear  gains  to  give  che  desired  range  dependence,  as  portrayed  in  Fig.  3.6-1. 
The  present  study  provides  the  first  CADET  application  in  which  the  range  de¬ 
pendent  components  of  the  seeker  noise  have  been  modelled  as  nonlinear  func¬ 
tions  of  state  variables.  In  previous  models,  the  missile -target  separation  was 
taken  to  be  linearly  decreasing  to  zero  at  the  nominal  terminal  time,  t^. 

The  boresight  error  distortion  model  includes  both  aberration  caused 
by  the  passage  of  incident  radiation  through  a  protective  cover  prior  to  detec¬ 
tion,  and  boresight  error  limiting  due  to  a  signal  processing  nonlinear  charac¬ 
teristic  which  may  be  introduced  to  circumvent  the  null  and  spurious  sidelobe 
response  inherent  to  the  restricted  beamwidth  of  the  antenna  or  detector.  The 
aberration  is  generally  a  complicated  effect,  expressed  as  an  angular  error, 

0ab’  a<^e<* to  seelcer  l°ok  angle,  0jook,  the  latter  being  the  angle  between 
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R-11602 


WHITE 

NOISE 

INPUTS 


Figure  3. 6-1  Nonlinear  Seeker  Noise  Model 


the  line-of-sight  to  the  target  and  the  missile  body  axis  defined  in  Fig.  3.6-2. 
The  aberration  effect  is  treated  extensively  in  Ref.  3  (refer  also  to  Fig. 

A.  6-2),  and  thus  is  not  considered  in  detail  here.  A  description  of  the  model 
employed  for  the  nonlinear  receiver  characteristic  is  given  in  Section  A.  6.  A 
functional  diagram  of  the  boresight  error  model  is  portrayed  in  Fig.  3.6-3. 


R-11600 


y-AXIS 


MISSILE 

CENTER  OF  G?  ilY 


*Hgure  3. 6-2  Seeker  Head  Configuration 
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R-U601 


SEEKER 

NOISE 


MEASURED 

BORESIGHT 

ERROR 


Figure  3. 6-3  Boresight  Error  Distortion  Model 


The  seeker  track  loop,  shown  in  Fig.  3.6-4,  is  designed  to  main¬ 
tain  the  measured  boresight  error  near  zero.  A  control  torque,  Tc,  is 
generated  which  has  a  component  proportional  to  the  boresight  error,  and 
suitable  dampling  is  assured  by  means  of  rate  feedback  provided  by  a  rate 
gyro  mounted  on  the  seeker  head.  A  compensator  of  the  proportional  plus 
integral  form  is  included  in  the  track  loop  to  remove  the  effect  of  steeiy 
state  disturbance  torques.  The  characteristics  of  the  control  loop  dynamics 
are  discussed  in  detail  in  Sections  A.  6  and  A.  7. 


H'  116  85 


SEEKER 


Figure  3.6-4  Seeker  Track  Loop  Model 
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The  seeker  head  disturbance  torque  model,  depicted  in  Fig.  3. 6-5, 
shows  the  three  effects  under  consideration.  The  most  complex  nonlinear 
phenomenon  is  an  effective  torque  due  to  seeker  mass  imbalance,  T  ,  with 
three  components;  one  proportional  to  missile  body  angular  acceleration, 
another  to  the  square  of  the  missile  body  angular  rate,  and  the  third  to  the 
missile  lateral  acceleration,  as  derived  in  Section  A.  6.  Another  nonlinear¬ 
ity  represents  an  applied  torque  Tr,  due  to  restoring  torques  acting  on  the 
seeker  head  (caused  either  by  wiring  harnesses  or  by  restraining  springs 
provided  to  prevent  large  seeker  head  angular  excursions).  Generally,  this 
is  a  "hard  spring"  effect,  i.e. ,  the  restoring  torque  will  be  negligible  for 
small  values  of  seeker  head  angle,  but  as  0^  increases,  the  restoring  torque 
will  increase  more  rapidly,  in  a  nonlinear  manner,  as  depicted  in  Fig. 
3.6-5.  This  effect  is  represented  here  by  a  power  law  nonlinearity,  viz. 


FROM 

AUTOttOT 


Figure  3. 6-5  Disturbance  Torque  Model 
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Tr  ’  W  ’  kl 


ev 

lim 


sign(6h) 


(3.6-1) 


where  k  is  an  integer.  A  final  nonlinearity  is  introduced  to  describe  the  fric¬ 
tion  torque,  T^;  in  particular,  we  consider  the  "dry"  or  coulomb  type  of  fric¬ 
tion  which  produces  a  counter  torque  that  is  constant  in  value,  having  the  sign 
of  the  seeker  angular  rate, 

Tf  =  VV  =  k2  Sign  (3.6-2) 


Coulomb  friction  can  cause  limit  cycles  to  appear  in  the  seeker  track  loop, 
producing  a  deterioration  in  the  missile  guidance  system  performance. 


3. 7  SUMMARY 

In  Fig.  3.7-1  the  complete  missile -target  intercept  model  is  portrayed 
with  all  of  the  subsystems  described  in  the  previous  sections  appropriately  in¬ 
terconnected.  All  of  the  system  variables  are  depicted  except  angle  of  attack, 
control  fin  deflection  and  the  seeker  compensation  state,  which  are  encompassed 
in  the  linear  dynamics  represented  by  the  transfer  functions  g1(s),  g2(s),  and 
(1  +kQ/s).  The  functions  f^Xj),  fg(x2)  and  fg(c)  represent  effects  due  to  non¬ 
linear  restoring  torques  acting  on  the  seeker  head,  nonlinear  friction  in  the 
seeker  gimbal,  and  the  receiver  characteristic,  respectively.  The  first  two 
are  described  in  Eqs.  (3.6-1)  and  (3.6-2).  The  last  function  is  the  standard 
limiter  operating  on  measured  boresight  error,  as  illustrated  in  Fig.  3.6-4. 

All  other  effects  are  describee  by  the  specific  functional  relations  shown. 
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4.  CADET  AND  MONTE  CARLO  STUDIES 

OF  NONLINEAR  EFFECTS 


4.1  CADET-MONTE  CARLO  COMPARISONS  FOR  THE  BASIC  SYSTEM 

In  extending  the  validity  of  CADET  as  a  missile  system  analysis 
tool,  we  first  consider  the  basic  system  model,  which  includes  the  follow¬ 
ing  nonlinear  effects: 

•  Range -dependent  seeker  noise  sources  (Fig.  3.6-1) 

•  Proportional  navigation  law  with  secant  compensation 
(Fig.  3.5-1) 

•  Acceleration  command  limiting  (Fig.  3.5-1) 

•  Missile-target  kinematics  (Section  3.1) 

•  Inverse  tangent  calculation  of  LOS  angle 

The  basic  system  model  incorporates  nine  distinct  nonlinearities,  listed  in 
Section  A.  8,  as  compared  with  previous  investigations  (Refs.  1  to  3)  where 
the  same  effects  were  approximated  by  only  three  nonlinearities  and  several 
linear  time -varying  gains. 

In  all  studies,  the  three  components  of  the  wide-band  seeker  noise 
(100  rad/sec  bandwidths)  are  always  specified  such  that  they  have  an  equal 
effective  rms  level  at  about  the  midpoint  of  the  terminal  homing  phase.  To 
be  more  precise,  the  spectral  densities  of  the  three  white  noise  inputs  w^  to 
w0  in  Fig.  3.6-1  are  chosen  to  achieve  rms  values  of  the  noise  states  x^, 
x^g  and  x^g  that  satisfy 
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q,  =  CT2  x  104  - CT. ,  =  0.7071  a  rad-ft 

1  o  14  o 

q0  =  a2  x  10"12  - ►  a1K  =  0.7071  x  10"8a  rad/ft 

z  o  Id  o 


0.7071  x  10~4ct  rad 
o 


(4.1-1) 


where  we  have  used  the  relation  indicated  in  Eq.  (3.3-2).  The  dimension¬ 
less  parameter  o  is  designated  the  seeker  noise  factor.  The  three  com¬ 
ponents  of  the  range  dependent  seeker  noise, 

ns  =  x14/r(t)  +x15r(t)  +  x16 

thus  have  rms  levels  given  approximately  by  cr^/ r(t),  r(t)  and  o^g, 

respectively,  which  are  equal  for  r(t)  =  10, 000  ft.  With  the  nominal  param¬ 
eter  values  indicated  in  Table  4.1-1,  r(t)  achieves  this  value  at  about 
t  =  2. 5  sec  in  a  terminal  homing  phase  of  about  six  seconds  duration. 

For  all  of  the  studies  that  follow,  the  fundamental  measure  of  sys¬ 
tem  performance  is  rms  miss  distance;  this  is  defined  to  be  the  rms  cross¬ 
range  missile -target  separation  at  the  mean  terminal  time,  tf.  The  mean 
terminal  time,  defined  as  the  instant  when  the  mean  down-range  separation, 
m  ,  goes  through  zero,  is  a  variable  in  this  study;  its  value  without  random 
disturbances  or  initial  conditions  is  6  sec,  while  the  presence  of  noise  sources 
or  a  nonzero  rms  target  maneuver  level  typically  increases  tf  by  a  few  hun¬ 
dredths  of  a  second.  This  effect  is  a  direct  result  of  relaxing  the  assumption 
of  constant  closing  velocity  that  was  made  in  earlier  work  (Refs.  1  and  2). 

In  every  case  considered  in  this  chapter,  one  or  several  of  the 
system  parameters  are  allowed  to  vary  from  the  nominal,  and  the  effects 
of  this  change  on  the  missile  guidance  system  performance  are  analyzed  via 
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TABLE  4.1-1 

NOMINAL  SYSTEM  PARAMETERS 
(Refer  to  Fig.  3.7-1) 


i  Seeker  head  mass,  m 
Seeker  head  moment  of  inertia  about  eg,  Iq 
Seeker  head  eg  offset,  tq 
Angle  of  seeker  head  eg  offset,  8g 
SEEKER  Distance,  missile  eg  to  seeker  gimbal,  rj 

PARAMETERS  Servo  gain,  ks 

Rate  gyro  gain,  kg 
Integral  compensator  gain,  kg 
Track  loop  time  conslant,  Tj 
Restoring  torque,  fj(x^) 

Gimbal  fricUan,  f2<X2) 


GUIDANCE 

LAW 

PARAMETERS 


AUTOPILOT/ 

AIRFRAME 

TRANSFER 

FUNCTIONS 


VEHICLE 

VELOCITIES 


RANDOM 

INPUT 

DANDWIDTHS 


NOISE  STATE 
RMS  LEVELS 


STATE 

VARIABLE 

INITIAL 

CONDITION 

STATISTICS 


Noise  filter  time  constant, 
Navigation  ratio,  n' 

Acceleration  command  limit,  amax 


¥S>  =  XT~2 


CgS  fCgS+Cj 


S  4CgS  ICgS+Cj 


d„8+d. 

c  (S)  .  _ - 2__1 - 

''  S  HCgS  -ICgS+C^ 


Missile  velocity  magnitude,  vm 
Target  velocity  magnitude,  v^ 


Seeker  noise  source  bandwidlhs,  coj,  o^,  Ug 
Target  maneuver  bandwidth, 

Range  rale  uncertainty  bandwidth,  (04 


Range  rate  uncertainty  level,  a ^ 
Target  acceleration  level,  o13 
Seeker  noise  factor  (Eq.(4.1-l)),c0 


Mean  heading  error,  mg 
RMS  heading  error,  erg 
Mean  down-range  separation,  m^j 
RMS  down-range  separation,  ojx 


0.15  oz-sec2/in 

0.1  in-oz-seefyrad 


40  In 

6  in-oz-sec/rad 


0.12  sec 
0  in-oz 
0  in-oz 
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275  sec*2 
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0.24  rad/ft -sec2 
O.C42  rad/ft-sec 
720  sec*3 
-0.865  sec*2 
-1.87  sec** 


3000  ft/sec 
1000  ft/sec 


100  rad/sec 
0.2  rad/sec 
100  rad/sec 


0  ft/sec 
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CADET.  A  few  of  the  resulting  performance  projections  are  compared  with 
monte  carlo  studies  to  provide  a  basis  for  assessing  the  accuracy  of  CADET 
in  capturing  the  effect  under  consideration.  Since  the  nominal  case  is  to  be 
used  as  the  starting  point  for  all  subsequent  studies  of  nonlinear  effects  and 
parameter  variations,  it  takes  on  a  significance  that  makes  it  advisable  to 
perform  a  monte  carlo  verification  using  a  large  number  of  trials. 

One  comparison  of  the  results  of  CADET  and  monte  carlo  analysis 
for  the  nominal  case  is  shown  in  Fig.  4.1-1.  A  second  case  with  a  reduced 
acceleration  command  limit,  amax  =  150  ft/sec^,  is  also  portrayed;  the 
latter  curve  demonstrates  the  significant  deterioration  in  the  missile  guid¬ 
ance  system  performance  that  results  from  reducing  the  missile  maneuver¬ 
ability.  This  nonlinear  effect  has  been  treated  in  previous  studies; 

Fig.  4.1-2,  taken  from  Ref .  2,  indicates  the  typical  variation  in  rms  miss 
distance  with  changing  missile  acceleration  limit. 


on  Basic  System  Performance 
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Figure  4.1-2  Acceleration  Command  Limit  Study  (Ref.  2) 


Returning  to  the  nominal  case  (a„ov  =  250  ft/sec^  in  Fig.  4.1-1), 
we  note  that  the  CADET  result  (the  smooth  curve,  exhibiting  an  rms  miss 
distance  of  22.9  feet  at  6  sec)  lies  within,  or  close  to,  the  95%  confidence 
band  of  the  corresponding  monte  carlo  data  throughout  the  engagement.  The 
width  of  the  confidence  band  provides  a  measure  of  the  theoretical  relia¬ 
bility  of  the  monte  carlo  data,  based  on  the  number  of  trials  performed  and 
an  estimate  of  X  (kurtosis;  refer  to  Section  2.1). 


The  fact  that  the  monte  carlo  estimate  of  the  rms  lateral  separa¬ 
tion,  a  ,  is  a  random  variable  that  converges  to  the  true  value  slowly  as 

y 

the  number  of  trials  increases  is  depicted  in  Fig.  4. 1-3,  where  we  show 
the  value  of  versus  the  number  of  trials  performed,  q.  The  case  shown 
in  Fig.  4.  l-3a  corresponds  to  the  estimation  of  rms  lateral  separation  two 
seconds  before  intercept,  when  the  probability  density  function  is  nearly 
gaussian  (X  ^  3).  The  accuracy  of  the  result  obtained  by  performing  500 
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Figure  4.  i-3 


(b)  t  =  6  sec 

Comparison  of  CADET  and  Monte  Carlo  rms 
Lateral  Separation  for  the  Nominal  Case 


monte  carlo  trials  is  quantified  by 

Pr  |^137. 7  ft  *  ay(4)  *  156.0  ft j  =  0.95 
Near  the  terminal  time,  however,  the  density  of  y  is  found  to  be  significantly 
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nongaussian;  on  the  basis  of  the  monte  carlo  study,  it  is  estimated  that  the 
kurtosis  is  approximately  15.  Thus  the  estimate  of  the  rms  miss  distance, 
0y(6),  necessarily  is  less  certain  than  when  y  is  gaussian.  The  large  value 
of  X  is  due  to  nonlinearity  in  the  system  which  gives  rise  to  a  small  but 
significant  probability  that  the  miss  distance  may  be  very  large  (see  also 
Section  5.3);  whenever  a  trial  is  performed  which  results  in  such  a  patho¬ 
logical  miss  there  is  a  large  ’'transient”  in  the  estimated  rms  miss  dis¬ 
tance.  This  phenomenon  is  observed  between  trials  140  and  180  and  between 
trials  420  and  460  in  Fig.  4.1-3b. 

The  agreement  of  CADET  and  the  monte  carlo  analysis  is  excellent 
as  long  as  y  is  reasonably  gaussian;  at  the  end  of  the  engagement,  CADET 
appears  to  have  underestimated  the  rms  miss  distance  by  ten  to  fifteen 
percent.  This  study  provides  a  graphic  demonstration  of  the  fact  that  sev¬ 
eral  hundred  or  even  a  few  thousand  monte  carlo  trials  are  required  to  ob¬ 
tain  an  estimate  of  rms  miss  distance  that  is  more  accurate  than  the  result 
given  by  CADET.  Taking  the  five  sets  of  100  trials  separately*,  the  value 
of  CTy(6)  varies  from  19.72  to  35.88  ft;  thus  the  use  of  only  100  trials  is  in¬ 
adequate  for  determining  the  absolute  accuracy  of  CADET.  On  the  basis  of 
this  quite  typical  behavior,  w?  generally  consider  that  the  CADET  result  is 
well  verified  whenever  it  lies  close  to  or  within  the  95%  confidence  band. 


The  present  study  of  the  nominal  case  also  shows  that  CADET  can  yield 
a  very  accurate  analysis  of  the  mean  value  of  the  down-range  separation, 


The  first  100  trials  shown  in  Fig.  4.  l-3b  yielded  cry(6)  -&i  -  19.72  ft,  the 
second  set  (trials  101  to  200)  resulted  in  cy(6)  ^0%-  32.08ft,  and  subse¬ 
quent  sets  yielded  <73  =  22.25  ft,  04  =  25.67  ft,  and  =  35.88  ft.  To 
aggregate  the  statistics  for  various  sets  of  trials,  it  is  necessary  to  aver¬ 
age  the  variances: 
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mx,  when  the  terminal  time  is  not  deterministic.  In  Fig.  4.1-4,  it  is 
demonstrated  that  the  CADET  ana  lysis  is  well  verified  by  the  monte  carlo 
results.  The  mean  terminal  time  given  by  CADET,  6.026  sec,  is  in  almost 
exact  agreement  with  the  monte  carlo  figure. 


Figure  4.1-4  Mean  Down-Range  Missile-Target  Separation 
Near  the  Mean  Terminal  Time 


The  remaining  two  components  of  the  missile-target  separation  are 

of  secondary  importance.  The  mean  lateral  separation,  my,  computed  by 

CADET  and  the  monte  carlo  method,  is  shown  in  Fig.  4.1 -5a.  We  note  that 

CADET  gives  the  exact  result  for  the  mean  of  y,  viz.  m  s  0  throughout  the 

** 
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(□)  CADET  AND  MONTE  CARLO  STUDY  Of  MEAN  LATERAL  SEPARATION 


(b)  CADET  AND  MONTE  CARLO  STUDY  OF  RMS  DOWN-RANGE  SEPARATION 

Figure  4 . 3.-5  Mean  Lateral  Separation  and  rms  Down 
Range  Separation  in  the  Nominal  Case 
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engagement,  *  while  the  monte  carlo  results  converge  to  a  zero  mean  quite  slowly. 
However,  CADET  incorrectly  indicates  that  x has  no  random  component  (a  =  0)  as 
shown  in  Fig.  4.  l-5b,  which  is  not  true,*  it  is  observed,  however,  that  is  negli¬ 
gible  except  in  the  last  few  hundredths^!  a  second  of  the  terminal  homing 
phase .  The  lack  of  a  random  component  of  x  in  the  CADET  analysis  is  due 
to  the  even  symmetry  of  the  nonlinearities  cos  (0a)  and  cos  (0^)  that  consti¬ 
tute  the  two  contributions  to  x;  any  even  function  with  a  zero  mean  input  has 
a  random  input  describing  function  gain  that  is  identically  equal  to  zero. 

Figure  4.1-6  indicates  the  sensitivity  of  the  miss  distance  to  changes 
in  the  noise  level  for  the  basic  system  configuration.  The  more  stringent 


Figure  4.1-6  Effect  of  Seeker  Noise  Level  on 
the  Basic  System  Performance 


* 

This  can  be  demonstrated  by  noting  that  there  are  no  mean  components  in 
the  initial  conditions  except  for  x(0)  =  24,000  ft,  and  no  bias  components  in 
the  random  disturbances,  so  the  basic  guidance  loop  nonlinearities  have 
odd  symmetry  and  all  variables  except  range  have  zero  means  for  all  time. 
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acceleration  command  limit,  a  =  150  ft/sec2,  was  retained  in  this  inves- 
tigation.  For  a  seeker  noise  factor  of  or  =  90,  a  significant  deterioration 
in  the  rms  miss  distance  is  observed,  and  the  CADET  analysis  is  well  veri¬ 
fied  by  the  monte  carlo  results. 

The  above  investigation  of  the  capability  of  CADET  to  provide  an 
accurate  measure  of  tactical  missile  performance  verifies  the  applicability 
of  CADET  to  the  basic  system  model,  thus  paving  the  way  for  studies  of 
other  effects. 


4.2  SEEKER  MASS  IMBALANCE 

In  this  study,  it  is  assumed  that  the  seeker  head  center  of  gravity 
(eg)  is  offset  from  the  gimbal  pivot  point  by  a  distance  rQ,  and  located  at  an 
angle  0q  with  respect  to  the  antenna  centerline  (refer  to  Fig.  3.6-2).  This 
eg  offset  results  in  an  effective  disturbance  torque  of  significant  complexity. 
Two  types  of  nonlinear ities  that  have  not  previously  been  studied  are  re¬ 
quired  to  model  the  effect  of  mass  imbalance.  These  are  the  product  of  one 
variable  with  a  trigonometric  function  of  another  variable  (vj_  cos  V2),  and 
the  square  of  a  variable  times  a  trigonometric  function  of  a  second  variable 
(vj  sin  v2). 

The  sensitivity  shown  in  Fig.  4.2-1  corresponds  to  the  situation 
where  the  seeker  head  eg  is  directly  in  front  of  the  gimbal  point  —  0Q  =  0  deg. 
Without  seeker  dynamic  compensation,  described  in  Section  3.6,  a  pro¬ 
nounced  deterioration  in  the  missile  guidance  system  performance  is  noted 
with  very  small  values  of  eg  offset.  Even  for  an  offset  of  one  thousandth 
of  an  inch,  Fig.  4.2-1  indicates  that  the  rms  miss  distance  is  nearly  twice  the 
nominal  value.  The  introduction  of  proportional  plus  integral  compensation 
to  remove  the  effect  of  steady  state  disturbance  torques  (Section  A. 7) 
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Figure  4.2-1  Effect  of  Geeker  Mass  Imbalance 
on  Guidance  Accuracy 


markedly  improves  the  guidance  system  performance.  In  both  the  uncom¬ 
pensated  and  compensated  cases,  the  monte  carlo  results*  are  in  good 
agreement  with  the  CADET  analysis,  even  when  extremely  large  miss  dis¬ 
tances  are  incurred.  Thus  CADET  gives  an  accurate  assessment  of  the 
effect  of  seeker  mass  imbalance . 


$ 

The  confidence  bands  plotted  in  Fig.  4.2-1  and  subsequent  figures  are 
based  on  the  gaussian  assumption  (X  s:  3).  In  the  previous  section  it  is 
demonstrated  that  the  actual  confidence  bands  are  probably  larger;  it 
is  not  possible  to  obtain  a  useful  estimate  of  the  kurtosis  with  fewer  than 
one  hundred  monte  carlo  trials,  however,  so  no  attempt  is  made  to 
consider  this  factor . 
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The  effect  of  the  angular  coordinate  of  the  eg  location,  0Q,  was 
studied  for  a  single  value  of  eg  offset  (0.03  in)  and  with  the  seeker  assumed 
to  be  compensated.  The  CADET  analysis  indicates  only  a  slight  difference 
in  rms  miss  distance  for  0q  =  180  deg  compared  with  the  case  0Q  =  0  deg, 
as  shown  in  Fig.  4.2-2.  With  the  eg  offset  normal  to  the  antenna  centerline 
(0q  =  90  deg),  the  guidance  system  terminal  performance  deterioration  was 
lessened  to  a  great  extent.  We  observe  that  the  most  significant  difference 
between  these  three  cases  occurs  in  mid-flight;  the  fact  that  the  maximum 
rms  lateral  separation  is  much  smaller  for  0q  =  0  deg  than  for  0Q  =  180  deg 
can  be  attributed  to  an  effective  change  in  the  seeker  dynamics. 


TIME,  t  (wt) 

Figure  4.2-2  Effect  of  Seeker  eg  Offset  Angie 
with  Compensated  Seeker 
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4.3  SEEKER  GIMBAL  DRY  FRICTION 

The  occurrence  of  dry  or  Coulomb  friction  in  the  seeker  gimbal 
introduces  a  disturbance  torque  of  the  form 

T.  =  k2  sign  (6h) 

where  &  is  the  seeker  angular  rate  measured  relative  to  the  missile  air- 
h 

frame.  Results  from  a  study  of  this  effect  are  presented  in  Fig.  4.3-1. 


SEEKER  DRV  FRICTION  COEFFICIENT,  kj  (in.-oz) 

Figure  4.3-1  Effect  of  Seeker  Gimbal  Dry 
Friction  on  Guidance  Accuracy 

First,  several  values  of  the  dry  friction  coefficient,  kg,  were 
selected  and  CADET  analysis  was  performed  to  obtain  a  prediction  of  the 
rms  miss  distance  that  would  result  without  seeker  dynamic  compensation. 
When  proportional  plus  integral  compensation  was  incorporated  in  the 
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seeker  track  loop,  a  significant  improvement  was  observed  in  the 
performance  projections  provided  by  CADET.  It  can  be  seen  in  Fig.  4. 3~1 
that  for  a  dry  friction  coefficient  of  0 . 8,  the  miss  distance  is  reduced  by  a 
factor  of  6  when  compensation  is  implemented. 

Three  monte  carlo  studies  were  undertaken  to  check  the  accuracy 
of  CADET;  two  for  the  compensated  case  and  one  to  verify  the  large  rms 
miss  distance  given  by  CADET  for  the  uncompensated  seeker.  In  every  case, 
the  CADET  result  was  near  or  within  the  95%  confidence  band  of  the  monte 
carlo  estimate,  which  verifies  the  ability  of  CADET  to  capture  the  effect  of 
dry  friction  quite  adequately  even  when  the  guidance  system  performance 
(as  measured  by  rms  miss  distance)  was  degraded  by  an  order  of  magni¬ 
tude  from  the  nominal  case. 

It  is  well  known  that  the  dry  friction  relay-type  characteristic  can 
lend  to  limit  cycling  in  dynamic  systems.  This  phenomenon  can  be  explained, 
using  the  standard  sinusoidal  input  describing  function  theory,  by  the  fact  that 
the  effective  gain  of  this  discontinuous  nonlinearity  is  extremely  large  for 
small  signals  (Ref.  8).  Thus  a  point  of  considerable  interest  is  whether  or 
not  the  missile  performance  degradation  observed  in  Fig.  4.3-1  is  due  to 
limit  cycles  in  the  seeker  track  loop.  A  single  monte  carlo  simulation  was 
performed  for  both  the  compensated  and  uncompensated  seeker,  with  dry  fric¬ 
tion  coefficients  of  1. 6  id  0. 8  respectively,  and  limit  cycles  were  clearly 
observed,  as  shown  for  the  compensated  case  in  Fig.  4.  3-2.  It  is  evident  that 
CADET  can  provide  a  reasonably  accurate  assessment  of  the  missile  guid¬ 
ance  system  performance  under  such  circumstances. 
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Figure  4.3-2  Existence  of  Limit  Cycles  due  to  Seeker  Gimbal 
Dry  Friction  with  Compensated  Seeker 


4.4  NONLINEAR  SEEKER  HEAD  RESTORING  TORQUES 

The  classical  "hard  spring"  type  of  nonlinearity  —  usually  modeled 
as  a  power  law  characteristic, 


0r 

lim 


sign  (0h) 


where  k  is  an  integer  greater  than  one  --  provides  an  example  of  another 
class  of  nonlinear  effects,  where  the  output  is  small  for  small  seeker  angle 
deflections  (0^  «  0^m)  but  increases  more  rapidly  than  a  linear  charac¬ 
teristic  when  0^  exceeds  Such  a  restoring  torque  might  occur  as  the 
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result  of  a  seeker  head  spring  restraint  system  designed  to  prevent  large 
values  of  seeker  head  angle  to  preclude  the  seeker  head  reaching  its  stops, 
or  might  occur  naturally  due  to  wiring  harnesses  or  other  flexible  linkage 
to  the  missile  body.  The  case  investigated,  k  =  11,  corresponds  to  a  negli¬ 
gible  spring  effect  for  small  0 h  and  a  very  rapid  increase  in  restoring  tor¬ 
que  as  0^  approaches  9lim-  In  this  study,  we  confined  our  attention  to  the 
compensated  seeker  case,  and  chose  0^.  to  be  10  deg. 

The  variation  of  rms  miss  distance  as  a  function  of  the  nonlinear 
spring  coefficient,  kp  is  shown  in  Fig.  4.4-1.  The  CADET  analysis  indi¬ 
cates  that  the  miss  distance  increases  very  abruptly  with  k^  for  k^  <0.2; 
the  monte  carlo  simulations  carried  out  for  two  cases  well  verified  the 
performance  projections  given  by  CADET. 


ft*TT0t1 


Figure  4.4-1  Effect  of  Nonlinear  Restoring  Torques 
on  Guidance  System  Performance 
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An  interesting  aspect  of  the  above  comparison  is  that  the  monte 
carlo  analysis  required  a  smaller  integration  step  size  in  simulating  the 
system  response  than  was  required  in  the  CADET  propagation  of  the  sta¬ 
tistics.  *  The  CADET  program  ran  correctly  with  an  integration  step  size 
of  0.01  sec,  while  the  monte  carlo  program  required  0.005  sec  to  avoid 
overflow  due  to  numerical  instability.  It  is  reasonable  to  conjecture  that 
CADET  may  be  more  immune  to  numerical  errors  inherent  to  the  inte¬ 
gration  techniques  for  solving  differential  equations  on  a  digital  computer, 
since  quasi-linearization  tends  to  "smooth"  nonlinearities  —  e.g.,  refer 
to  Fig.  5.1-2;  thus  in  cases  where  the  nonlinearity  exhibits  large  vari¬ 
ations  in  slope,  CADET  will  tend  to  have  an  increased  computational  ad¬ 
vantage  with  respect  to  the  monte  carlo  method. 

4. 5  RECEIVER  CHARACTERISTIC  AND  RANGE  RATE  UNCERTAINTY 

The  effect  of  the  receiver/signal  processing  nonlinearity  for  the 
case  where  the  boresight  error  limiter  (Fig.  3.6-3)  saturates  at 
elim  =  ^-25  deg  (representative  of  narrow-beamwidth  monopulse  radar 
systems  (Ref.  12)  or  infrared  seekers)  was  considered.  The  effectiveness 
of  the  seeker  was  reduced  by  increasing  the  track  loop  time  constant  from 
0.12  to  0.32  sec  (which  increases  the  boresight  error  rms  level  by  the 
same  factor,  2.6,  and  thus  increases  the  effect  of  limiting).  For  the 
above  parameter  values,  Fig.  4.5-1  shows  that  CADET  predicted  a  deteri¬ 
oration  in  rms  miss  distance  from  22.9  ft  (in  the  nominal  case)  to  28. 8  ft. 
The  monte  carlo  method  provided  a  good  verification  of  this  result;  for 
60  trials,  the  nominal  rms  miss  distance  was  predicted  to  be  20.9  ft  and 
the  effect  of  choosing  =  0.25  deg  and  =  0. 32  sec  was  to  increase 

* 

Both  programs  used  the  same  integration  technique  —  the  fourth-order 
Runge-Kutta  method. 
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Figure  4.5-1  Effect  of  Boresight  Error  Saturation  on  Guidance  Accuracy 

this  miss  to  27.2  ft,  where  the  same  sets  of  initial  conditions  and  noise 
inputs  were  used  for  the  two  cases . 

The  effect  of  range  rate  uncertainty  in  the  guidance  law  (refer  to 
Fig.  3. 5-1)  was  assessed  for  the  case  of  a  fixed  (bias)  uncertainty.  It  was 
assumed  that  the  guidance  package  range  rate  information  was  in  error  by 
a  fixed  value  of  500  ft/sec.  The  monte  carlo  result,  shown  in  Fig.  4.5-2, 
is  that  the  rms  miss  distance  is  decreased  slightly  by  this  error;  CADET 
analysis  predicts  that  the  effect  is  virtually  insignificant  (Fig.  4.5-3).  Both 
methods  agree  that  the  rms  lateral  separation  is  reduced  from  the  nominal 
case  in  mid-flight  for  e  ^  =  +500  ft/sec;  this  is  to  be  expected,  since  the 


THE  ANALYTIC  SCIENCES  CORPORATION 


R-12013 


Figure  4.5-2  Effect  of  Bias  Range  Rate  Uncertainty 
on  Guidance  System  Performance, 

20  Trial  Monte  Carlo  Analysis 

error  modeled  above  is  equivalent  to  assuming  that  the  guidance  package 
range  rate  value  is  approximately  4500  ft/sec  rather  than  about  4000  ft/sec 
in  the  nominal  case.  *  This  corresponds  to  a  range  rate  error  of  +  12.  5%, 
which  is  equivalent  in  effect  to  the  same  percentage  increase  in  the  navi¬ 
gation  ratio,  n'.  Previous  studies  (Ref.  2)  have  shown  that  the  guidance 
system  performance  is  quite  insensitive  to  moderate  changes  in  n7;  how¬ 
ever,  a  higher  value  of  n7  results  in  a  somewhat  faster  guidance  loop  which 
should  lead  to  the  reduction  in  rms  ’ateral  separation  in  mid -flight. 

Closing  velocity  is  actually  a  random  variable,  as  indicated  in  Fig.  3.5-1. 
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TIME,  i  (sec) 

Figure  4.5-3  Effect  of  Bias  Range  Rate  Uncertainty  on 
Guidance  System  Performance  via  CADET 

4.6  CADET-MONTE  CARLO  COMPARISON  WITH  ALL 
NONLINEAR  EFFECTS 

As  a  final  study  of  the  efficacy  of  CADET  in  providing  an  accurate 
statistical  analysis  of  the  performance  of  a  complex,  highly  nonlinear  tac¬ 
tical  missile  guidance  system  model,  we  consider  the  effect  of  exercising 
all  of  the  nonlinearities  in  combination.  Values  of  the  parameters  are 
shown  in  Table  4.6-1.  They  were  chosen  such  that  each  nonlinearity  alone 
had  led  to  a  significant  deterioration  of  system  performance  from  nominal 
in  the  studies  treated  in  Sections  4.2  to  4 . 5 . 

The  results  of  the  CADET  and  monte  carlo  statistical  analysis  are 
shown  in  Fig.  4.6-1.  The  CADET  value  of  rms  lateral  separation  is  well 
verified  by  the  30 -trial  monte  carlo  study. 
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TABLE  4.6-1 

PARAMETER  VALUES  THAT  EXERCISE  ALL  NONLINEARITIES* 


Effect 

Parameters 

Values 

Seeker  Mass  Imbalance 

r0’  0O 

0.03  in,  0  deg 

Nonlinear  Restoring  Torque 

kl'  "<  6l,m 

0.05  in-oz,  11,  10  deg 

Nonlinear  Gimbal  Friction 

k2 

0.8  in-oz 

Boresight  Error  Limiter 

clim 

0.25  deg 

Seeker  Compensation  Gain 

k0 

20.0 

All  Nonlinearities  Exercised 

% 

Only  deviations  from  the  nominal  case,  Table  4.1-1,  are  noted. 
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4. 7  COMPARISON  OF  CADET  AND  MONTE  CARLO  EFFICIENCY 

One  of  the  main  arguments  advanced  for  the  use  of  CADET  in  ob¬ 
taining  projections  of  tactical  missile  guidance  system  performance  is  the 
significant  reduction  in  computer  central  processing  unit  (CPU)  time  achieved 
by  CADET  when  compared  with  the  monte  carlo  method  of  statistical  analysis. 

In  making  this  evaluation,  two  issues  must  be  addressed:  the  number  of  monte 
carlo  trials  that  must  be  performed  in  order  to  obtain  comparably  accurate 
results,  and  the  practical  limitation  imposed  by  computer  costs.  From  the 
standpoint  of  accuracy,  any  such  number  is  somewhat  arbitrary,  because  the 
error  mechanisms  of  CADET  and  the  monte  carlo  method  are  essentially  dis¬ 
similar,  and  of  necessity  the  judgement  is  based  on  a  limited  body  of  experi¬ 
ence.  Referring  to  Fig.  4.1-3b,  we  note  that  in  a  situation  where  the  statistics 
are  quite  nongaussian  the  CADET  computation  of  rms  lateral  separation  appears 
to  be  at  least  as  accurate  as  the  value  estimated  with  400  monte  carlo  trials,  in 
the  sense  that  the  95%  confidence  band  for  400  trials  brackets  the  CADET  result. 
Where  the  statistics  are  more  nearly  gaussian,  e.g.  as  in  Fig.  4.1-3a,  it 
would  appear  that  one  thousand  or  more  trials  may  be  required  in  order  to 
achieve  comparable  accuracy.  On  the  other  hand,  a  pragmatic  evaluation  of 
the  efficiency  of  CADET  should  take  into  account  the  fact  that  most  monte  carlo 
studies  must  be  limited  in  scope  by  computer  budget  constraints.  A  reasonable 
upper  bound  is  256  trials  since,  in  the  gaussian  case,  this  results  in  95%  con¬ 
fidence  that  an  accuracy  of  10%  can  be  achieved  (Section  2.1);  for  high  order 
systems,  even  this  number  of  trials  may  require  too  much  computer  time. 

For  the  present,  we  will  thus  compare  the  relative  efficiency  of  the  monte 
carlo  and  CADET  approaches  on  the  basis  of  256  trials,  recognizing  that  the 
estimated  rms  values  of  the  system  variables  obtained  for  this  number  of 
monte  carlo  experiments  may  be  less  accurate  than  the  CADET  results. 

In  the  present  study,  the  savings  in  computer  CPU  time  is  quite  sig¬ 
nificant,  even  though  the  system  is  of  considerably  higher  order  (n  =  17)  and 
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has  more  nonlinearitios  than  in  previous  investigations  that  compared  CADET 
and  the  monte  carlo  technique.  Both  of  these  factors  tend  to  reduce  the  relative 
efficiency  of  CADET.  Monte  carlo  simulation  requires  only  the  integration  of 
an  n-vector  differential  equation  (q  times),  while  CADET  involves  the  propagation 
of  an  n-vector  and  annxn  symmetric  matrix  --  a  total  of  n  (n  +  3)/2  elements. 
Thus  the  computational  burden  for  CADET  can  increase  as  fast  as  n  /2  while 
the  CPU  time  for  monte  carlo  analysis  only  varies  as  n,  demonstrating  that  an  in¬ 
crease  in  the  number  of  states  may  reduce  the  advantage  of  CADET  in  efficiency. 
This  factor  can  be  mitigated  where  there  is  little  cross-coupling  in  the  system; 
in  the  quasi-linear  system  model,  this  corresponds  to  having  few  non- zero 
elements  (Nf  being  sparse).  In  many  practical  problems,  is  sparse  and  a 
considerable  increase  in  the  computational  efficiency  of  CADET  can  be  realized 
by  the  application  of  techniques  which  circumvent  multiplications  involving  zero 
elements,  thus  streamlining  the  evaluation  of  P  (Eq.  (2.2-5)). 

The  number  of  nonlinearities  may  also  increase  the  computation  time 
required  by  CADET,  since  the  calculation  of  a  random  input  describing 
function  in  CADET  generally  requires  more  logical  and  numerical  opera¬ 
tions  than  evaluating  the  corresponding  nonlinear  function  in  the  monte 
carlo  program  (refer  to  Appendix  B).  The  present  study  was  exceptional 
in  having  nearly  as  many  nonlinearities  as  states;  more  typical  applica¬ 
tions  of  CADET  would  focus  on  a  few  principal  nonlinear  effects,  leading 
to  a  greater  reduction  in  computational  burden  per  performance  evalua¬ 
tion  in  comparison  to  a  monte  carlo  analysis. 

Using  the  same  integration  method  in  performing  the  monte  carlo 
ensemble  of  simulations  as  was  used  in  propagating  the  system  mean  vector 
and  covariance  matrix  via  CADET,  and  assuming  that  the  same  integration 
step  size  can  be  used  in  p'  on  procedure,  it  has  been  possible  to  perform 
10  CADET  sensitivity  studies  at  the  same  computational  expense  needed 
for  one  accurate  monte  carlo  study.  Since  we  have  seen  that  it  is  possible 
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that  the  monte  carlo  approach  may  require  a  reduced  integration  step  size 
to  avoid  failure  of  the  numerical  integration  technique  (Section  4.4),  the 
ratio  may  be  even  higher.  For  the  case  studied  in  Section  4.4,  the  CADET 
analysis  required  only  4.4  minutes  of  CPU  time,  in  contrast  to  about  92 
minutes  for  a  comparable  256  trial  monte  carlo  study  with  one-half  the  step 
size.  A  more  typical  application  of  CADET  can  exhibit  even  greater  ad¬ 
vantages;  a  simpler  but  still  realistic  missile -target  intercept  problem  has 
been  treated  in  which  a  ratio  of  30  CADET  analyses  to  one  monte  carlo 
study  was  achieved  (Ref.  2). 

The  studies  discussed  in  this  chapter  have  treated  the  missile- 
target  intercept  problem,  represented  by  a  nonlinear  model  of  considerable 
complexity  —  with  between  9  and  16  nonlinearities.  The  results  presented 
demonstrate  that  CADET  has  been  quite  successful  in  capturing  the  significant 
sensitivities  of  rms  miss  distance  to  the  nonlinear  phenomena  in  question. 
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5 .  SENSITIVITY  STUDIES  AND  COMPUTATIONAL  ISSUES 


5.1  EFFECTS  OF  PROBABILITY  DENSITY  FUNCTIONS  ON 
RANDOM  INPUT  DESCRIBING  FUNCTIONS 


An  important  issue  that  must  be  investigated  in  order  to  assess  the 
potential  success  of  CADET  in  providing  accurate  performance  projections  for 
tactical  missile  guidance  systems  is  the  effect  of  the  assumption  that  the 
state  variables  are  jointly  normal  on  the  calculation  of  random  input  describ¬ 
ing  functions  (ridf's).  The  gaussian  hypothesis  is  the  only  basic  approxima¬ 
tion  made  in  the  application  of  CADET,  so  any  inaccuracy  in  the  statistical 
analysis  obtained  via  CADET  is  due  to  deviation  of  the  actual  joint  probability 
density  function  (pdf)  from  normality. 


In  this  section,  the  sensitivity  of  CADET  to  changes  in  the  pdf  of  the 
nonlinearity  input  is  investigated  by  comparing  the  ridf’s  corresponding  to 
selected  nonlinearities  often  found  in  missile  guidance  system  models,  com¬ 
puted  for  a  variety  of  density  functions.  Three  nonlinearities  are  chosen; 
these  are  the  limiter,  the  sinusoidal  operator,  ana  a  power  law  nonlinearity. 
Seven  probability  density  functions  with  quite  different  functional  forms  are 
considered.  Four  of  these  are  taken  from  Table  2.1-1,  viz,,,  the  exponential, 
gaussian,  triangular,  and  uniform  distributions.  Three  additional  densities 
are  special  cases  of  the  sum  of  two  symmetrical  triangular  functions,  gen¬ 
erally  defined  by 


P(x) 


1 

2  A 


0 


!  |  X  |  -  XQ|  ^  A 

(5.1-1) 

|  |  x  |  -  xJ>A 
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wh:  Jh  has  a  zero  mean,  a  variance  given  by 


2  2  1.2 
a  =  x  +  ^  A 
o  o 


and  a  ratio  of  fourth  moment  to  variance  squared  of 


,4  "4 


4.22  1.4 

o  o  15 

T  1a2 . 2 — FTl 

x  +tA  x  +««A 
o  3  o  36 


(5.1-2) 


(5.1-3) 


The  three  cases  of  Eq.  (5.1-1)  chosen  for  the  present  study  correspond  to 
A  =  ~xq,  XQ  and  2xq;  the  associated  ndf*s  are  portrayed  in  Fig.  5.1-1.  Note 
that  two  of  these  densities  are  bimodal;  i.e.,  they  have  ;o  peaks. 
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wh:  Jh  has  a  zero  mean,  a  variance  given  by 


2  2  1  a2 
a  =  x  +  ~-  A 
O  b 


and  a  ratio  of  fourth  moment  to  variance  squared  of 


4  *2  2  1.4 

x  +  A  x  +  vc-A 
o  o  15 

“i  1.2"  2 — T~i 

x  +  tA  x^  +  oc  A 
o  3  o  36 


(5.1-2) 


(5.1-3) 


The  three  cases  of  Eq.  (5.1-1)  chosen  for  the  present  study  correspond  to 
A  =  -Xq,  x0  and  2xq;  the  associated  pdf's  are  portrayed  in  Fig.  5.1-1.  Note 
that  two  of  these  densities  are  bimodal;  i.e.,  they  have  /o  peaks. 
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in  Eq.  (5.1-4)  ha va  even  symmetry,  we  note  that  skew  densities  can  be  dis¬ 
regarded  with  no  loss  in  generality.  For  a  skew  density  p  (x)  we  can  define 

s 

its  even  part  by 


Pey(x)  =  |  (PsW  +  Ps("x)) 

Since  the  three  nonlinearities  considered  are  symmetric  (odd)  and  the 
mean  values  of  their  inputs  are  zero,  only  the  even  part  of  the  pdf  contributes 
to  the  describing  function  calculation. 

Limiter  -  The  ideal  limiter, 

Ix  ,  \x\z  6 

(5.1-5) 

6  sign  (x)  ,  |xl  •>  6 

is  probably  the  most  common  piece -wise  linear  function  used  to  model  non¬ 
linear  phenomena;  here  it  represents  several  saturation  effects  in  the 
missile  guidance  system.  In  Fig.  5.1-2,  we  portray  the  various  describing 
function  gains  for  this  nonlinearity,  corresponding  to  the  pdf’s  defined  in 
Eq.  (5.1-4),  as  functions  of  the  ratio  of  the  input  rms  level,  o,  to  the  satura¬ 
tion  point,  6.*  As  would  be  expected,  all  seven  quasi-linear  gains  capture  the 
fact  that  the  effective  gain  starts  to  decrease  from  unity  whenever  a  significant 
portion  of  the  assumed  input  pdf  lies  beyond  the  saturation  point,  i.e.,  when¬ 
ever  there  is  a  significant  probability  that  |xl  is  greater  than  6.  As  has  been 
pointed  out  previously,  this  effect  is  the  key  to  the  success  of  quasi-lineariza¬ 
tion  techniques  in  reflecting  nonlinear  system  behavior  that  is  beyond  the  scope 
of  small-signal  (Taylor  series)  linearization. 


* 

The  derivations  of  these  and  all  subsequent  ridf's  are  given  in  Appendix  B. 
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R  »l**? 


Figure  5.1-2  Random  Input  Describing  Function  Sensitivity  for 
the  Limiter 


It  is  interesting  to  observe  that  the  relative  position  of  the  curves 
in  Fig.  5.1-2  exhibits  a  monotonic  relation  to  the  value  of  X.  The  greater 
the  difference  between  X  for  a  particular  pdf  and  the  value  for  the  gaussian 
case  (X  =  3),  the  greater  the  difference  between  that  density  function's  ridf 
curve  and  the  curve  for  a  gaussian  distribution.  This  behavior  holds  in  all 
the  cases  considered  here,  and  is  indicative  of  the  fact  that  the  value  of  Xis 
one  quantitative  measure  of  how  "close"  the  density  function  is  to  being 
gaussian. 

The  variation  of  the  ridf's  with  X  is  about  at  its  maximum  (on  a  per¬ 
centage  basis)  for  the  case  cr  *  26.  This  is  shown  in  Fig.  5.1-3;  we  note 
that  the  ridf  decreases  13%  as  X  increases  from  3  to  6,  and  it  increases  28% 
as  X  decreases  from  3  to  1,16. 
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R-11971 


RATIO  OF  FOURTH  CENTRAL  MOMENT  TO 
THE  VARIANCE  SQUARED,  X 

Figure  5.1-3  Random  Input  Describing  Function  Sensitivity  to 
X  for  the  Limiter,  a  =  2  6 


Power  Law  -  A  similar  study  was  performed  for  a  power-law 
characteristic, 

f(x)  =  x2  sign  (x)  (5.1-6) 

This  type  of  nonlinearity  is  often  used  to  model  effects  such  as  the  "hard 
spring"  characteristic,  as  treated  in  the  study  of  nonlinear  restoring  torques 
acting  -n  the  seeker  head.  For  the  power  law,  the  ridf's  calculated  for  the 
same  density  functions  considered  previously  are  shown  to  be  of  the  form 

nr  =  M.cr  (5.1-7) 

i 
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in  Appendix  B,  where  a  is  the  input  rms  level  and  are  coefficients  deter¬ 
mined  by  the  input  pdf's,  p.fx).  Thus  the  describing  function  gain  for  f(x) 
increases  linearly  with  the  input  rms  level,  in  direct  contrast  to  the  small- 
signal  linear  gain  which  is  identically  equal  to  zero,  as  shown  in  Fig.  5.1-4. 

It  is  again  observed  that  there  is  a  monotonic  relation  between  X^  and  the  ridf 
curves,  depicted  by  plotting  u.  versus  X^  in  Fig.  5.1-5.  In  this  case,  an  in¬ 
crease  in  X  leads  to  an  increase  in  the  describing  function  gain,  which  is  con¬ 
trary  to  the  behavior  shown  for  the  limiter.  This  is  a  result  of  the  fact  that 
the  power  law  output  increases  more  rapidly  with  increasing  input  than  a 
linear  characteristic,  whereas  the  opposite  is  true  for  the  limiter.  For  the 
power  law  nonlinearity,  the  ridf  sensitivity  is  independent  of  a,  i.e.,  the 

ratio  of  ridf’s  calculated  for  p.(x)  and  p.(x)  is  simply  u./m..  For  f(x)  in 

^  1  ^  3 

Eq.  (5.1-6),  the  gain  nr  varies  from  +33%  for  the  exponentially  distributed 
case,  to  -34%  for  the  pdf  p^(x),  compared  to  the  gaussian  input  ridf,  which 
shows  that  this  nonlinearity  is  somewhat  more  sensitive  to  variations  in  X 
than  the  limiter . 

Sinusoidal  Operator  -  The  third  nonlinearity  considered  in  these 
sensitivity  studies  is  the  sinusoidal  operator, 

f(x)  =  sin  x 

which  represents  the  resolution  of  the  missile  and  target  velocity  vectors 
into  orthogonal  components  in  the  missile -target  intercept  mode].  A  potential 
source  of  difficulty  with  this  function  is  that  the  nonlinearity  output  periodically 
changes  sign  with  increasing  or  decreasing  values  of  its  input.*  This  leads 
to  quasi- linear  gains  that,  for  large  values  of  rms  input,  o,  may  even  differ 
in  sign  for  different  input  pdf’s.  This  problem  is  not  unique  to  CADET;  in 
many  modeling  and  simulation  studies,  care  must  be  exercised  - 1  .he  input 
to  a  sinusoidal  operator  can  exceed  +  90  deg  (+  ^  rad),  s*-  sense 


♦This  is  not  a  problem  in  the  present  study. 
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Figure  5.1-4 


Random  Input  Describing  Function  Sensitivity  for 
the  Power  Law  Nonlinearity 


EXPONENTIAL 

'GAUSSIAN 
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RATIO  OF  FOURTH  CENTRAL  MOMENT  TO 
THE  VARIANCE  SOUARED.  X 


Figure  5.1-5 


Random  Input  Describing  Function  Sensitivity 
to  \  for  the  Power  Law  Nonlinearity 
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the  "gain"  can  become  negative  in  this  situation.  Bearing  this  in  mind,  we 
have  calculated  the  random  input  describing  functions  for  values  of  cr  as  large 
as  3  rad  to  indicate  where  such  effects  become  important,  as  shown  in 
Fig.  5.1-6. 


INPUT  STANDARD  DEVIATION,  0  (r*d) 


Figure  5.1-6  Random.  Input  Describing  Function  Sensitivity  for 
the  Sinusoidal  Operator 


The  quasi -linear  gains  for  a  <  |  rad  show  some  similarity  to  those 
obtained  for  the  limiter;  this  is  a  reasonable  mode  of  behavior,  since  the  sine 
function  shows  a  definite  saturation  effect  over  the  range  |x|  s  |  rad.  As 
expected,  the  ridf's  are  inversely  related  to  X  for  a<  ^  rad,  i.e.,  as  X  in¬ 
creases,  nr  decreases,  as  shown  in  Fig.  5.1-7  for  a  =  1.0  rad.  However, 
as  the  input  rms  level  approaches  3  rad,  the  describing  functions  for  all  of 
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the  pdf's  except  p^(x)  and  p£(x)  become  negative,  and  the  monotonic  relation¬ 
ship  between  X  and  nr  appears  to  be  lost. 

The  preceding  studies  show  that  in  all  three  cases  the  sensitivity  of 
random  input  describing  function  calculations  to  variations  in  input  proba¬ 
bility  density  function  is  slight  for  small  values  of  input  rms  level;  as  0 
approaches  zero,  the  quasi-linear  gains  approach  unity  for  the  limiter  and 
sinusoidal  operator,  and  zero  for  the  power  law  nonlinearity.  These  limit¬ 
ing  cases  are  the  same  values  of  gain  that  would  be  obtained  by  the  tradi¬ 
tional  small-signal  linearization  approach  --  viz. ,  by  replacing  f(x)  with  a 
linear  gain  equal  to  the  slope  of  the  function  at  the  origin  (Section  2.2).  As 
a  general  result,  it  has  been  shown  (Ref.  8)  that  quasi -linearization  sub¬ 
sumes  small-signal  linearization,  i.e. ,  for  small  signals  the  two  are 
equivalent.  This,  in  turn,  proves  that  CADET  provides  nearly  exact  sta¬ 
tistical  analyses  when  the  random  variables  have  a  small  rms  value  in 
relation  to  the  system  nonlinearities,  i.e.,  when  most  of  each  nonlinearity 
input  probability  density  function  lies  in  the  linear  region  of  its  nonlinearity. 
As  the  rms  levels  of  system  variables  increase  so  that  the  nonlinearities 
are  being  exercised  significantly,  the  describing  function  sensitivity  to  the 
input  pdf  can  be  appreciable;  then  it  must  be  ascertained  how  sensitive  the 
system  performance  is  to  variations  in  gain  at  each  point  in  the  system 
model  where  a  nonlinearity  occurs.  No  general  answer  can  be  given  to 
this  question;  the  verification  of  CADET  for  particular  applications  must 
be  accomplished  by  direct  comparison  with  monte  carlo  results,  as  has 
been  done  in  Chapter  4  for  the  missile  homing  guidance  system. 
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5.2  RANDOM  INPUT  DESCRIBING  FUNCTIONS  NOT  EXISTING 
IN  CLOSED  FORM  UNDER  THE  GAUSSIAN  ASSUMPTION 

For  certain  nonlinearities,  the  question  of  sensitivity  has  a  strong 
bearing  on  computational  issues.  These  are  cases  for  which  the  random  in¬ 
put  describing  functions  cannot  be  obtained  in  closed  form  under  the  assump¬ 
tion  that  the  states  are  jointly  normal.  An  example  of  interest  in  the  missile- 
target  intercept  problem  is  the  nonlinearity 

f 2  2 

r  =yjx  +  y 


which  defines  the  missile -to-target  range  in  terms  of  the  down-range  and 
cross-range  components,  x  and  y,  respectively.  In  order  to  simplify  the  dis¬ 
cussion,  we  assume  that  y  does  not  have  a  mean  component,  and  x  has  a 
negligible  random  component  except  at  the  very  end  of  the  engagement  (refer 
to  Section  4.1).  This  approximation  is  valid  in  the  missile -target  intercept 
situations  investigated  here. 

In  this  case,  there  are  only  two  ridf’s,  nmx  and  nry,  needed  for  a 
quasi-linear  representation  of  the  range .  Thus  we  must  evaluate 


dy 


dy 


(5.2-1) 


Under  the  assumption  that  y  is  a  gaussian  random  variable,  the  second  of  these 
integrals  can  be  evaluated  analytically;  however,  the  first,  which  is  of  the  form 
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(5.2-2) 


cannot  generally  be  solved  in  closed  form  unless  m  =  0,  viz., 

A 


e[^P]  =  E[|yl3  =71  ay  (5.2-3) 


For  the  more  general  case  given  by  Eq.  (5.2-2)  with  m  -  0  it  is  desirable 

A 

to  use  some  approximate  technique  to  obtain  a  closed  form  expression  for 
nm  that  is  convenient  for  use  in  a  CADET  analysis.*  The  two  principal 
approaches  that  we  consider  are  approximation  of  the  nonlinearity  (e.g., 
by  series  expansion),  and  substitution  of  approximate  pdf's. 

A  Taylor  series  expansion  of  a  function  of  a  random  variable,  f(y), 
about  the  mean  of  that  variable,  here  taken  to  be  zero,  results  in 


£(y)  =  f(0)  + 


/df_\  1  ( d2f\ 

W/y=0y+2\dy2iy=O 


(5.2-4) 


We  desire  to  determine  the  expected  value  of  the  function,  which  is  given  by 


E[f(y,]  =  £<0)+H$)y=0ffy2+H$)-0E[y3]+‘"(5‘2'5) 

% - 

While  n  in  Eq.  (5. 2-2)  can  be  calculated  by  numerical  integration,  a  less 
x 

time  consuming  approach  is  desired  for  repeated  evaluation  in  a  CADET 
analysis. 
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where  use  is  made  of  the  fact  that  E[y]  is  zero  to  eliminate  the  second  term  in 
Eq.  (5.2-4);  all  other  odd-moment  terms  (E  [y^]  etc.)  are  also  zero  for 
symmetric  pdf’s.  Truncating  the  series  given  in  Eq.  (5.2-5)  at  the  second 
term,  we  obtain 


E[f(y)]  er  f(0)  + 


'A\ 

dy2/y=0<T 


(5.2-6) 


which  is  an  approximation  suggested  in  Ref.  9.  We  note  that  this  result  is 
independent  of  the  particular  density  function  of  y.  If  more  terms  are  de¬ 
sired,  the  higher-order  moments  can  be  evaluated  using  a  specified  pdf.  If 
y  is  gaussian,  all  odd  moments  are  zero  and  even  moments  are  given  by  (Ref  .9) 


"2k  =  E 


r  2k*| 

y 


(1)  (3)  (5) 


(2k  -  1)  a2k  ,  k  =  1,  2, 


which  leads  to  the  full  expansion 


„!■*/  m  t/n\  l/d2f\  2  (1)  (3)/ d  i\  4 
E[f(y)3  =  f(0)+o(— 2’)  CT  — T  +. 

2  \dy  /  y=0  \dy  /  y-0 


(5.2-7) 


The  use  of  the  first  term  alone  in  Eq.  (5.2-7)  corresponds  to  small 
signal  linearization;  taking  two  terms  as  indicated  in  Eq.  (5.2-6)  results  in  a 
quasi-linear  gain  that  is  often  useful.  In  the  present  case,  however,  this 
approach  is  effective  only  in  situations  where  m  is  considerably  larger  in 

A  — 

magnitude  than  a  ,  due  to  the  singularities  of  the  derivatives  of  Jm  2  +  y2 
y  v  a 

at  the  origin  (m  =  0).  To  demonstrate  this  difficulty,  we  write  the  series 
expansion  for  the  nonlinearity  under  consideration  (Ref.  10), 
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f(y)  =  j 


m  +  y 
x  3 


|  m  | 
1  x' 


1  + 


Hki-m 


from  which  we  obtain 


(5.2-8) 


B 


m 


1  + 


2 .  i  /°l\ 

2  U  )  8  (m  l 


4 


sign(mx)  (5.2-9) 


as  an  approximate  describing  function  to  represent  the  mean  component  of  the 
range.  For  mx  considerably  larger  than  the  first  few  terms  of  this  ex¬ 
pansion  yield  acceptable  accuracy.*  Since  m  approaches  zero  in  the  missile- 
target  intercept  problem,  however,  this  result  is  generally  not  suitable. 

Another  approximation  to  the  nonlinearity  leads  to  a  series  that  is 
useful  for  small  values  of  m  : 

X 

^mx2  +  y2  -  lyl  +  !mxje“^y//mx'  ,  ImJ  <  I y |  (5.2-10) 


The  right-  and  left -hand-sides  of  Eq.  (5.2-10)  have  the  same  first  three  terms 
when  their  power  series  expansions  are  compared.  Substituting  this  approxima¬ 
tion  into  Eq.  (5. 2-2)  leads  to  an  integral  that  can  be  evaluated  in  terms  of  the 
complementary  error  function,  denoted  as  erfc  (Ref.  11), 


We  note  that  the  expansion  indicated  in  Eq.  (5.2-9)  never  converges  formally, 
i.e.,  for  any  value  of  ay/mx,  no  matter  how  small,  the  series  will  eventually 
diverge  as  more  terms  are  evaluated.  This  is  a  standard  property  of  asymp¬ 
totic  expansions  which  are  useful  only  when  truncated  after  a  finite  number 
of  terms. 
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which  in  turn  can  be  represented  as  an  asymptotic  expansion  (Ref.  6)  valid 
for  or  »  | m  l: 

y  x 


For  m  =  0,  the  expansion  correctly  yields  the  result  given  in  Eq.  (5.2-3), 

y.  j 

and  for  |m  |  <  ■=■  a  ,  the  result  is  quite  accurate.*  However,  the  series  is 
x  y 

divergent  for  larger  values  of  |m  | .  This  series  is  of  questionable  utility 

A 

in  the  missile -target  intercept  problem,  since  |m  |»cr  at  the  beginning  of 

X  y 

the  terminal  homing  phase . 

The  second  method  for  approximately  evaluating  the  first  integral 
in  Eq.  (5.2-1)  is  the  substitution  of  a  nongaussian  pdf  for  which  the  integral 
can  be  obtained  in  closed  form.  As  in  previous  sensitivity  studies  (Section  5.1), 
the  best  result  has  been  obtained  using  the  triangular  pdf.  Substituting  this 
distribution  into  Eq.  (5.2-1)  leads  to  an  integral  that  is  evaluated  in  closed 
form  to  be 


(5.2-12) 


*Eq.  (5.2-11)  is  also  not  formally  convergent. 
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where  the  auxiliary  variable  v  is  given  by 


v 


(5.2-13) 


The  accuracy  erf  Eq.  (5.2-12)  is  quite  good,  especially  when  compared  with 
the  poor  approximations  given  by  the  series  expansion  in  Eqs.  (5.2-9)  and 
(5.2-11)  when  |mx|  is  nearly  equal  to  ffy .  The  error  between  Eq.  (5.2-12) 
and  the  exact  result  specified  in  Eq.  (5.2-2)  is  less  than  3%,  which  is  ade¬ 
quate  for  most  applications. 


In  passing,  we  note  a  third  alternative  that  might  be  considered  in 
applications  where  less  accuracy  is  required;  viz. 


sign  (mx) 


(5.2-14) 


The  simplicity  of  this  result  —  and  the  fact  that  the  approximation  is  better 

than  either  series  expansion  (Eqs.  (5.2-9)  and  (5.2-11))  for  |m  |  ^rj  — 

x  y 

makes  it  attractive,  despite  the  error  of  25%  at  m  =0. 

All  of  the  approximate  solutions  for  nmx  in  Eq.  (5.2-2)  discussed 
above  are  compared  in  Fig.  5.2-1.  Generally,  we  note  that  se  es 
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I 

i 

! 


approximations  for  this  nonlinearity  are  accurate  only  in  applications  where 
i  mxj » cjy  or  1  raj  <<  cry  for  the  entire  CADET  simulation;  such  conditions 
are  never  satisfied  in  studying  the  expected  value  of  the  range  in  the  missile 
target  intercept  problem.  The  use  of  an  approximate  (nongaussian)  pdf  — 
in  particular,  the  triangular  density  —  yields  the  best  overall  accuracy. 


Ill 

o 

2 

< 

cc 

111 

X 

cc 

o 

u. 

2 

O 

H 

O 

2 

D 


s 

o 

Q 
2 
< 
C C 

o 

U1 

N 

Zj 

< 

2 

tc 

o 

2 


R-11973 


Figure  5.2-1  Comparison  of  Approximations  for  the  Expected 
Value  of  the  Range 


It  should  be  mentioned  that  the  conclusion  that  the  series  expansion 
techniques  are  not  very  accurate  for  computing  the  ridf  in  the  case  treated 
above  should  not  be  taken  as  universally  true.  When  series  approximations 
for  an  ridf  can  be  obtained  which  are  accurate  over  the  entire  range  of  the 
input  statistics,  they  will  generally  yield  good  results. 
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Ultimately,  the  issue  of  sensitivity  depends  upon  how  much  the 
CADET  analysis  results  —  e.g.,  the  computed  rms  miss  distance  —  are 
affected  by  errors  or  variations  in  the  describing  function  gains.  It  is  con¬ 
ceivable  that  the  miss  distance  will  not  be  significantly  affected  by  large 
changes  in  a  particular  ridf,  depending  upon  how  the  ridf  enters  into  the 
system  model.  In  the  missile  guidance  system  model  developed  in  Appendix  A 
the  only  incidence  of  range  dependence  is  in  the  seeker  noise  module.  It  is 
intuitively  clear  that  the  effect  of  seeker  noise  on  the  statistical  performance 
of  the  guidance  system  will  not  be  particularly  sensitive  to  moderate  inaccur¬ 
acies  in  the  quasi-linear  representation  of  range,  particularly  if  these  in¬ 
accuracies  occur  only  during  the  last  few  hundredths  of  a  second  of  the  engage 
ment.  This  is  verified  by  comparing  CADET  results  obtained  for  the  nominal 
case  (Section  4.1)  using  the  two  range  ridf  approximations  labeled  ’’triangular 
density”  and  ”Eq.  (5.2-9),  2 Terms"  in  Fig.  5.2-1.  The  difference  in  the 
CADET  computed  rms  miss  statistics  for  these  cases  was  negligible;  in  fact, 
all  state  variables  had  means  and  standard  deviations  that  differed  by  less 
than  one  percent  at  terminal  time . 

On  the  basis  of  the  insensitivity  of  the  seeker  noise  module  to  varia¬ 
tions  in  the  ridf  representation  of  range  dependent  noise  sources,  we  recom¬ 
mend  the  use  of  a  simpler  ridf  than  that  obtained  from  the  triangular  density 
approximation  (Eq,  (5.2-12)).  Rather  than  using  the  two-term  series  expan¬ 
sion  considered  in  the  above  comparison,  however,  the  application  of  Eq. 
(5.2-14)  should  be  an  even  more  attractive  alternative,  due  to  the  fact  that 
the  latter  result  is  more  accurate  for  small  values  of  m  .  In  particular, 

Eq.  (5.2-14)  does  not  become  infinite  as  m  goes  to  zero. 

A 


5.3  HISTOGRAM  STUDIES  OF  NONGAUSSIAN  SYSTEM  VARIABLES 

In  considering  the  impact  of  the  assumption  that  the  state  variables 
are  jointly  normal  on  the  accuracy  of  performance  projections  provided  by 
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CADET,  we  have  performed  a.  histogram  study  of  a  esse  in  which  there  is  a 
significant  difference  between  CADET  and  monte  carlo  results.  Such  a  study 
provides  some  insight  into  how  nongaussian  the  state  variable  pdfs  must  be 
in  order  for  CADET  results  to  be  inaccurate. 

The  system  under  consideration  corresponds  to  the  nominal  case 
given  in  Table  4.1-1,  except  for  a  single  large  initial  condition  statistic  — 
an  rms  initial  heading  error  of  10  deg  —  chosen  to  obtain  a  significant  dis¬ 
crepancy  between  CADET  and  monte  carlo  results;  this  somewhat  unrealistic 
condition  leads  to  a  very  large  rms  miss  distance. 

The  time  history  of  the  rms  lateral  separation,  o  ,  over  the  duration 

J 

of  the  engagement  is  presented  in  Fig.  5.3-1.  The  monte  carlo  confidence 
limits  shown  are  calculated  on  the  preliminary  assumption  that  y  is  nearly 
gaussian,  i.e.,  that  X  sr  3  (Section  2.1).  The  CADET  and  monte  carlo  results 
agree  well  over  the  first  half  of  the  homing  phase.  However,  there  is  a 
marked  divergence  after  t  =  3  sec,  which  is  evidently  due  to  some  inadequacy 
in  the  assumption  that  the  states  are  jointly  normal. 

To  understand  the  reason  for  this  problem,  and  to  assess  how  much 
deviation  from  normality  is  required  for  a  discrepancy  of  this  sort  to  occur, 
we  portray  typical  monte  carlo  generated  histograms  for  a  number  of  the  sys¬ 
tem  states,  obtained  from  a  100-trial  ensemble  of  simulated  missile-target 
intercept  engagements.  First,  we  consider  the  target  aspect  angle,  0  ,  which 
is  theoretically  gaussian  over  the  entire  flight.  The  histograms  shown  in 
Fig.  5.3-2  correspond  to  the  midpoint  and  end  of  the  engagement.  We  note 
that  they  appear  to  be  quite  irregular;  this  is  due  to  the  fact  that  100  trials 
are  insufficient  for  obtaining  an  accurate  representation  of  the  probability 
density  function.  Thus,  for  example,  the  fact  that  there  appears  to  be  a  small 
"valley”  slightly  to  the  left  of  zero  at  t  =  6  sec  should  not  be  interpreted  to 
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mean  that  0  has  a  bimodal  distribution.  These  histograms  provide  the  basis 

A 

for  a  realistic  assessment  of  the  departure  of  other  histograms  in  this  study 
from  the  normal  case. 

Next,  we  consider  a  system  state  variable  that  is  distinctly  non- 
gaussian  for  this  trajectory,  the  missile  fin  deflection,  6 .  Because  of  the 
large  initial  rms  heading  error,  there  is  a  high  probability  that  the  missile 
initially  attempts  to  achieve  maximum  lateral  acceleration.  Thus  the  output 
of  the  acceleration  command  limiter,  ac  in  Fig.  3.5-1,  will  tend  to  be  bimodal, 
as  will  the  fin  deflection,  which  is  linearly  related  to  a  .  At  t  =  1  sec  a  highly 
bimodal  distribution  is  noted  in  Fig.  5.3-3.  As  the  terminal  homing  phase 
proceeds,  6  becomes  somewhat  less  bimodal,  but  at  no  time  does  it  appear 
to  be  gaussian. 

Other  variables  that  exhibit  bimodal  behavior  throughout  this  engage¬ 
ment  are  the  seeker  angle  and  angle  rate,  0^  and  0^;  missile  body  angle  and 
angle  rate,  and  0m;  and  angle  of  attack,  a.  We  observe  in  Fig.  5.3-4 
that  the  missile  lead  angle,  0^,  which  is  separated  from  6  by  a  significant 
amount  of  linear  dynamics  (refer  to  Fig.  3.7-1)  remains  nearly  gaussian 
throughout  the  engagement  —  thus  providing  an  example  of  the  "filter  hypothesis" 
that  the  presence  of  several  stages  of  low-pass  filtering  tends  to  make  a 
signal  more  nearly  gaussian,  which  forms  the  basis  for  the  gaussian  assump¬ 
tion  in  CADET. 

The  cross-range  missile -target  separation,  y,  exhibits  a  behavior 
that  more  directly  provides  an  explanation  of  the  failure  of  CADET  to  provide 
accurate  lateral  separation  statistics  over  the  second  half  of  the  terminal 
homing  phase.  For  the  first  half  of  the  engagement,  Fig.  5.3-5  indicates 
that  the  pdf’s  for  y  are  quite  nearly  gaussian.  However,  at  4  sec  we  note 


5-22 


NUMBER  OF  TRIALS 


THE  ANALYTIC  SCIENCES  CORPORATION 


THE  ANALYTIC  SCIENCES  CORPORATION 


fl 

n 

n 


TH^  ANALYTIC  SCIENCES  CORPORATICN 


.  J 


(a)  t=  4s«c 

Figure  5.3-6 
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Cb)  t  =  6  lie 


Histograms  for  Cross -Range  Separation  in  the 
Second  Half  of  the  Homing  Phase 


in  Fig.  5. 3 -6a  that  significant  "tails"  have  developed  in  the  histogram;  at 
6  sec  this  effect  is  more  pronounced  (Fig.  5.3-6b).  A  quantitative  com¬ 
parison  of  the  latter  histogram  with  the  gaussian  density  is  obtained  by  com¬ 
puting  the  ratio  of  the  fourth  central  moment  to  the  variance  squared;  this  is 
an  estimate  of  the  parameter  X  defined  in  Eq.  (5.1-3).  The  resulting  value 
of  19.5,  which  is  more  than  six  times  larger  than  the  value  for  a  gaussian 
density,  indicates  a  highly  nongaussian  distribution.  Two  conclusions  can  be 
drawn: 


•  The  monte  carlo  estimated  rms  miss  distance, 
CTy(tf)  =  1018  ft,  is  statistically  considerably  less 
accurate  than  implied  by  the  confidence  band  shown 
in  Fig.  5.3-1.  under  an  assumed  gaussian  distri¬ 
bution  for  y.  For  X=  19. 5,  the  95%  confidence  inter¬ 
val  limits  on  a  are 
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ff  =  750  ft,  a  =  2570  ft 

-y  y 


as  given  by  Eq.  (2.1-14). 

•  The  CADET  result  may  also  be  inaccurate,  since  the 
sensitivity  studies  in  Section  5.1  indicate  that  ridf 
calculations  based  on  the  gaussian  assumption  (X=  3) 
can  be  quite  different  from  those  based  on  a  pdf  with 
X  =  19.5. 

The  latter  point  indicates  the  reason  why  CADET  seems  to  provide  an  in¬ 
accurate  guidance  system  performance  projection  over  the  last  3  seconds 
of  flight. 


To  explain  some  of  the  above  observed  behavior,  an  analysis  of  the 
physical  significance  of  large  initial  heading  error  is  useful.  We  consider 
a  simplified  example,  with  the  aid  of  Fig.  5.3-7,  wherein  a  target  is  pro¬ 
ceeding  precisely  along  the  original  LOS  with  constant  velocity,  vt,  and 
not  maneuvering.  The  missile  has  been  launched  with  an  initial  heading 
error  of  0^  degrees,  with  constant  velocity  vector  magnitude,  vm«  An 
important  limitation  is  imposed  on  the  missile  by  assuming  that  lateral 
acceleration  is  constrained  by  jam 
command  limiting;  this  condition  makes  it  impossible  for  intercept  to 
occur  in  this  simple  scenario  if  |0^  |  exceeds  a  certain  value.  To  demon¬ 
strate  this  point,  we  assume  that  am  =  amax  over  the  entire  flight.  Then 
the  missile  will  follow  a  circular  trajectory  with  a  radius  and  angular  rate 
given  by  (Ref.  13) 


£  a  ,  as  a  result  of  acceleration 
max’ 


r  =  v  /a 
c  m  max 


0  =  a  /v 

c  max  m 


(5.3-1) 
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Figure  5.3-7  Intercept  with  a  Circular  Missile  Trajectory 

as  shown  in  Fig.  5. 3-7.  If  an  intercept  is  achieved  by  applying  the  maximum 
acceleration  over  the  entire  interval  0  £  t  *  t^,  the  equations  of  motion  are 

xt(t)  =  xto  -  vt‘ 


yt(t)  5  o 


xm(t)  =  rc[sine*0-sin<940-ec>] 
ym<t)  =  rc  [cos  (8i0  -  9C)  -  cos 


(5.3-2) 


(5.3-3) 


According  to  Eq.  (5.3-1), 


a  t 
max 


(5.3-4) 
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In  general,  in  order  for  intercept  to  be  possible, the  point  at  which  the 

circular  arc  followed  by  the  missile  intersects  the  x-axis  must  be  in  front  of 

the  target;  i.e.,  if  the  earliest  time  the  missile  trajectory  can  intersect  the  x- 

axis  is  t,,  it  muist  be  true  that 
r 

Xm(tf)  S  (5.3-5) 

to  permit  an  intercept.  The  values  of  0  .A  and  a  determine  whether  or  not 
Eq.  (5.3-5)  can  be  satisfied.  From  the  conditions  depicted  in  Fig.  5.3-7,  it  is 
clear  that  intercept  is  possible  only  for  i  I  ^  0£O>  where  and  the 
corresponding  terminal  time  t^  satisfy 


Xt0  ”  Vtfcf  =  Vm^  sin  9,o  -  sin  (ojto  ^ /amax 


cos^0=cos4l0--^-ij 


(5.3-6) 


These  transcendental  equations  have  the  solution 


®j40  =  14,4 


t<  =  6.048  sec 
f 


(5.3-7) 


for  the  parameters  used  in  generating  the  CADET  results  given  in  Fig.  5  3-1. 
For  9^q  greater  than  the  limit  0(,q  given  in  Eq.  (5.3-7),  the  deterministic 
miss  distance  (defined  in  this  context  to  be  the  absolute  value  of  ym(tf)) 
increases  very  rapidly,  as  shown  in  Fig.  5.3-8. 
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Figure  5.3-8  The  Variation  of  Miss  Distance  with  Initial 
Heading  Error  Magnitude 


From  the  above  analysis,  we  can  appreciate  why  the  rms  heading 
error  chosen  for  the  analysis  in  Fig.  5.3-1  leads  to  a  large  miss  distance. 
Since  0^  is  a  gaussian  random  variable  with  a  standard  deviation  of  10°, 
there  is  a  significant  probability  that  1 0^1  exceeds  the  maximum  value 


Prob  ^0^1  a  14.4  degj  =0.15 

and  the  impact  of  initial  heading  errors  greater  than  14.4  deg  in  magnitude  is 
dramatic  (Fig.  5.3-8).  In  fact,  we  can  apply  the  statistical  properties 
of  0&q  to  the  characteristic  shown  in  Fig.  5.3-8  to  derive  the  idealized 
rms  miss  distance  (valid  only  for  the  above  simplified  model  of  the  missile - 
target  intercept  problem),  defined  by 


00 

3y2f|E[y2(e10)]  =/  y2(et0)  P(el0)  dei0  (5.3-9) 

-  CO 
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The  resulting  value  of  is  736  feet.  W.s  can  also  perform  a  random  input 
describing  function  analysis  of  the  same  effect  using  Fig.  5.3-8  to  demon¬ 
state  why  CADET  has  given  an  estimated  rms  miss  distance  that  is  consider¬ 
ably  less  than  the  monte  carlo  result.  The  characteristic  shown  in  Fig.  5.3-8 
is  in  effect  a  linear  gain  with  a  dead  zone,  whose  random  input  describing 
function  given  by  Ref.  8  can  be  expressed  in  terms  of  the  complementary  error 
function,  erfc, 

n  =  309  erf c  =  48. 6  ft /deg  (5.3-10) 

r  \1Q^2/ 


Thus  the  rms  miss  distance  calculated  using  describing  function  theory  is 


=  488  feet 


(5.3-11) 


In  the  idealized  case  treated  above,  quasi-linearization  significantly  under¬ 
estimates  the  rms  miss  distance.  This  is  reflective  of  the  fact  that  the  ridf 
calculated  in  Eq.  (5.3-10)  is  based  on  the  minimization  of  mean  square 
approximation  error,  not  on  variance  matching. 

This  section  has  analyzed  a  missile-target  engagement  situation 
where  CADET  and  monte  carlo  results  are  significantly  different,  in  order  to 
provide  insight  as  to  potential  sources  of  error  in  the  CADET  results.  How¬ 
ever,  it  is  typically  found  that  parameter  values  must  be  chosen  which  generate 
large  rms  miss  distance  (more  than  100  ft)  in  order  to  obtain  noticeable  deter¬ 
ioration  in  CADET  accuracy.  In  those  cases  the  missile  dynamics  are  highly 
nonlinear  throughout  the  trajectory,  and  generally  do  not  correspond  to  condi¬ 
tions  under  which  tactical  missiles  are  designed  to  operate.  Furthermore, 
if  the  rms  miss  distance  is  large,  say  500  ft  or  more,  it  is  not  too  important 
that  CADET  makes  a  100  foot  error  if  the  missile  lethal  radius  is  50  feet,  since 
a  miss  distance  of  4Q0  ft  is  still  intolerable.  For  more  realistic  trajectories 
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that  yield  tolerable  ms  miss  distances,  CADET  results  appear  to  be  accu¬ 
rate  within  ±  10%,  based  upon  the  monte  carlo  simulations  performed  in  this 
and  other  studies. 

Finally,  we  observe  that  the  marked  deviation  of  the  density  of  the 
terminal  cross-range  missile -target  separation,  y,  from  the  gaussian  case 
is  a  potential  source  of  misinterpretation  of  the  value  of  a  .  The  cumulative 
distribution  of  ly  |  at  the  terminal  time  for  50  monte  carlo  trials  is  portrayed 
in  Fig.  5. 3-9;  the  actual  rms  value  of  y  for  this  data  is  0^  =  1018  ft.  On  the 
same  plot  we  show  cumulative  distributions  corresponding  to  |y  |  for  a 
gaussian  random  variable  y  having  cr^  =  1018  ft  (corresponding  to  the  monte 
carlo  result)  and  =  513  ft  (CADET).  In  terms  of  guidance  system  per¬ 
formance,  if  we  assume  that  it  is  desired  to  have  |y  |  less  than  50  feet  at 
the  terminal  time,  we  see  in  Fig.  5-  3-9  that  the  actual  distribution  indi¬ 
cates  an  effectiveness  of  60%,  i.e. ,  60%  of  the  engagements  satisfy  the  con¬ 
dition  |y  |  <  50  feet,  while  the  corresponding  values  of  effectiveness  for  the 
two  gaussian  distributions  are  8%  for  CADET  and  4%  for  monte  carlo  results, 
based  on  the  miss  distance  statistics  alone.  Thus  we  observe  that  the  rms 
miss  distance  is  itself  of  questionable  value  as  a  measure  of  guidance  sys¬ 
tem  performance  when  the  nongaussian  nature  of  y  leads  us  to  question  the 
accuracy  of  the  CADET  and  monte  carlo  results.  In  one  sense,  then,  the 
"failure"  of  CADET  in  this  situation  may  be  linked  to  the  "failure"  of  rms 
miss-distance  as  a  single -parameter  measure  of  system  performance. 
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6.  SUMMARY  AND  CONCLUSIONS 

6.1  SUMMARY 

The  principal  goal  of  this  study  is  to  extend  the  proven  capability 
of  the  Covariance  Analysis  DEscribing  Function  Technique  --  CADET  —  to 
provide  accurate  performance  projections  for  tactical  missile  guidance  sys 
tem  models  that  are  quite  realistic  —  i.e.,  that  incorporate  a  number  of  sig 
nificant  nonlinear  and  random  effects.  The  approach  used  to  achieve  this 
objective  has  entailed 

•  Verification  of  CADET  performance  projections  by 
the  use  of  selected  monte  carlo  performance  studies 

•  Investigation  of  the  sensitivity  of  CADET  analysis  to 
deviation  from  the  assumption  that  the  state  variables 
are  jointly  normal. 

In  this  investigation,  the  following  quite  diverse  effects  were  treated: 

Sources  of  Nonlinearity 

•  Guidance  law 

•  Acceleration  command  limiting 

•  Missile-target  intercept  geometry  (4  nonlinearities) 

•  Coordinate  transformation  (2  nonlinearities) 

•  Range -dependent  seeker  noise  sources  (2  nonlinearities) 

•  Receiver  /.signal  processing  characteristic 

•  Seeker  mass  imbalance  (3  nonlinearities) 


6-1 


THE  ANALYTIC  SCIENCES  CORPORATION 

•  Seeker  gimbal  Coulomb  friction 

•  Seeker  head  restoring  torques  (nonlinear  spring  effects) 

Random  Effects 

•  Tracking  sensor  noise  and  measurement  errors 

•  Range  rate  measurement  error 

•  Target  maneuvers 

•  Deviation  of  initial  conditions  from  nominal  values 

The  generation  of  missile  performance  statistics  via  CADET,  and 
the  subsequent  verification  of  selected  points  on  the  parameter  sensitivity 
curves  by  monte  carlo  simulation  techniques  is  treated  in  Chapter  4.  The 
basic  system  model  —  incorporating  the  first  five  nonlinear  effects  listed 
above,  which  are  fundamental  to  the  missile -target  intercept  problem  — 
was  studied  first,  then  each  of  the  remaining  nonlinearities  were  investi¬ 
gated  singly,  and  a  final  CADET -monte  carlo  comparison  was  made  with  all 
nonlinear  phenomena  present.  Three  aspects  of  the  sensitivity  problem  are 
considered  in  Chapter  5:  the  sensitivity  of  random  input  describing  function 
calculations  to  the  probability  density  function  of  the  nonlinearity  input,  calcu¬ 
lation  of  approximate  random  input  describing  functions  when  it  is  inconvenient 
to  use  the  result  for  the  gaussian  case,  and  a  histogram  study  of  a  situation 
where  system  variables  are  quite  highly  nongaussian. 

6.2  CONCLUSIONS 

The  investigation  described  in  this  report  has  indeed  shown  that 
CADET  is  an  accurate  and  efficient  tool  for  conducting  statistical  analyses 
of  the  performance  of  a  tactical  missile  system,  including  the  effects  of  a 
number  of  significant  nonlinear  and  random  phenomena.  The  conclusions 
drawn  from  the  study  can  be  summarized  as  follows: 
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•  CADET  has  the  demonstrated  ability  to  capture  the  effect  of 
each  of  the  16  nonlinearities  treated  in  Chapter  4  on  guidance 
system  performance.  In  all  cases  studied,  CADET  results 
are  close  to  or  within  the  95%  confidence  limits  of  the  monte 
carlo  analysis.  This  degree  of  agreement  was  generally  main¬ 
tained  even  in  the  numerous  instances  where  the  nonlinearities 
were  shown  to  have  a  marked  deleterious  effect  on  rms  miss 
distance. 

•  Despite  the  high  order  of  the  system  (17  state  variables)  and 
large  number  of  nonlinearities  (16),  CADET  shows  a  signifi¬ 
cant  computational  advantage  over  the  monte  carlo  method: 
Between  10  and  20  CADET  performance  projections  can  be 
obtained  for  the  same  amount  of  computer  time  required  by 
one  accurate  monte  carlo  study . 

•  It  is  sometimes  necessary  to  use  approximate  random  input 
describing  functions  due  to  difficulty  in  deriving  closed  form 
solutions.  Two  approaches  are  available,  based  on  series 
expansion  techniques  and  on  the  substitution  of  alternative 
density  functions  in  lieu  of  the  gaussian  density.  A  discus¬ 
sion  of  some  benefits  and  disadvantages  is  given  in  Chapter  5. 

•  Highly  nongaussian  system  variables  not  only  lead  to  inaccu¬ 
racy  in  the  CADET  analysis,  but  also  make  the  monte  carlo 
method  less  reliable  and  reduce  the  meaningfulness  of  the 
basic  measure  of  system  performance,  rms  miss  distance. 

•  The  value  of  the  parameter  X  (the  fourth  central  moment  of  a 
density  function  divided  by  the  variance  squared)  is  a  useful 
measure  of  the  departure  of  the  density  of  a  random  variable 
from  the  gaussian  case.  It  would  thus  be  valuable  to  estimate 
this  parameter  for  each  nonlinearity  input  in  the  monte  carlo 
analysis  to  help  in  appraising  the  accuracy  of  the  monte  carlo 
method  and  CADET. 


In  light  of  these  and  related  findings,  it  is  felt  that  confidence  in  the  applica¬ 
bility  of  CADET  to  perform  statistical  analyses  of  complex  nonlinear  missile 
guidance  systems  with  a  number  of  random  disturbances  has  been  significantly 
\  enhanced. 
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APPENDIX  A 

MISSILE -TARGET  INTERCEPT  MODEL 


A.  1  INTRODUCTION 

This  appendix  presents  the  derivation  of  the  mathematical  model  of 
the  general  missile -target  intercept  problem  used  for  verifying  the  capability 
of  CADET  to  provide  an  efficient,  accurate  statistical  analysis  of  guidance  sys¬ 
tem  performance.  The  e.  sential  guidelines  used  in  deriving  the  model  are 
realism,  completeness,  and  the  desire  to  include  enough  significant  nonlinear 
and  random  effects  to  establish  confidence  in  the  capability  of  CADET  to  treat 
complex  nonlinear  missile  systems  with  random  inputs. 

The  overall  interconnection  of  the  subsystems  which  comprise  the 
missile -target  intercept  model  is  indicated  in  Fig.  A.  1-1.  The  principal 
variables  are  shown  as  outputs  of  the  appropriate  blocks,  and  random  distur¬ 
bances  are  denoted  w..  Detailed  models  underlying  each  input-output  relation¬ 
ship  are  given  in  subsequent  sections  of  the  appendix.  While  the  basic  structure 
of  the  guidance  system  remains  essentially  unchanged  from  an  earlier  study 
described  in  Ref.  2,  there  are  numerous  extensions  and  refinements  incorpora¬ 
ted  in  the  present  investigation. 

The  final  model  is  a  state  vector  differential  equation  having  the  for¬ 
mulation 


x  =_f(x)  +  Gw(t)  (A.  1-1) 

where  w(t)  is  a  vector  of  gaussian  white  noise  inputs  which  represents  various 
random  effects,  the  matrix  G  determines  which  state  or  states  are  driven  by 
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MISSILE  *-HMO 


Figure  A.  1-1  Basic  System  Block  Diagram 


each  noise  input,  and  _f  (x)  embodies  the  linear  and  nonlinear  dynamic  relations 
within  the  system. 


A.  2  THE  MISSILE -TARGET  KINEMATICS  MODEL 

The  missile  target  engagement  treated  in  this  study  is  restricted  to  the 
terminal  homing  phase  in  a  planar  intercept  configuration.  Both  the  target  and 
the  missile  are  assumed  to  have  a  constant  magnitude  velocity  vector  which  may 
be  rotated  by  the  application  of  a  lateral  (normal)  acceleration. 
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An  initial  coordinate  system  is  defined  by  the  positions  of  the  missile 
and  target  at  the  initiation  of  the  terminal  homing  phase  (taken  to  occur  at 
t  =  0);  the  missile  is  at  the  origin  and  the  line-of -sight  (LOS)  to  the  target  de¬ 
fines  the  x-axis  at  t  =  0  (Fig.  A.  2-1).  The  origin  moves  with  the  missile, 
without  rotation,  so  that  x  and  y,  respectively,  provide  the  instantaneous  down 
range  and  cross  range  missile -target  separation.  Expressing  the  separation 
in  polar  coordinates,  the  relations 


V  2  2 

r  = !x  +  y 

8  =  tan"* (y/x) 


(A.  2-1) 


define  the  instantaneous  range  and  LOS  angle  of  the  target.  The  angles  0^ 
missile  lead  angle)  and  0  (target  aspect  angle)  specify  the  orientation  of  the 

cL 

missile  and  target  velocity  vectors;  by  convention,  6.  and  0  are  positive  in 

X/  a 

the  directions  defined  in  Figure  A.  2-1. 
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Figure  A.  2-1  Target -Mssile  Planar  Intercept  Geometry 
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In  terms  of  the  variables  defined  in  Fig.  A.  2-1,  the  differential 


equations 


•  a  •  1 

®  l  ~  x8  ~  v~  am 
m 


6  =  xQ  =  —  a, 
a  9  vt  t 


(A.  2-2) 


y=i10  =  vtsin<x9)-vmsin(x8) 


x  =*11=  -vt  cos  (x9)  -  vm  cos  (x8) 


express  the  dynamics  of  the  missile -target  separation  under  the  above  assump¬ 
tions,  and  identify  the  state  variables  (elements  of  x  in  Eq.  (A.  1-1))  xg  to  x^. 
The  LOS  angle  0  and  range  r  are  obtained  by  nonlinear  operations  on  x  and  y, 
given  in  Eq.  (A.  2-1)  and  shown  in  Fig.  A.  2-2. 
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Figure  A.  2-2  State  Vector  Formulation  of 
Missile-Target  Kinematics 
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A.  3  THE  TARGET  MODEL 

The  model  representing  the  target  behavior  remains  unchanged  from 
prior  studies  (Ref.  2).  As  mentioned  previously,  the  target  velocity  vector  is 
assumed  to  have  a  constant  magnitude,  and  a  direction  described  by  0  in 

2L 

Fig.  A.  2-1.  The  target  maneuver  model  used  in  this  study  represents  target 
lateral  acceleration  as  a  band -limited  gaussian  process  derived  from  a  gaussian 
white  noise  input  by  one  stage  of  low-pass  filtering.  In  differential  equation 
formulation,  we  have 


*13  =  '  “t  x13  +  w5  (A-  3_1> 

where  the  state  variable  is  the  target  lateral  acceleration,  a^;  the  equiva¬ 
lent  low -pass  filter  representation  is  depicted  in  Fig.  A.  3-1. 


* -1193V 
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Figure  A.  3-1  Target  Maneuver  Model 


The  target  maneuver  rms  level  ar^3  is  determined  by  the  spectral  den¬ 
sity  of  Wg  and  the  initial  condition  on  o-^; 


E[w5(t)  w5(t)]  =  qg  5(t  -  t) 
E[x13(0)2]  =  -gi 


(A.  3-2) 


t 


gives  us  an  rms  level  that  is  constant  throughout  the  engagement, 


CT13  s  V«5/2wt 


(A. 3-3) 
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It  is  interesting  to  note  that  the  autocorrelation  function  and  the  cor¬ 
responding  power  spectral  density  for  a  poisson  square  wave  --  i.  e  a  square 
wave  that  switches  between  i  a^ft/sec  with  random  poisson  distributed  switch¬ 
ing  times  having  an  average  of  o^/2  zero-crossings  per  second  —  are  identical 
to  those  of  the  above  gaussian  process,  although  the  associated  probability  den¬ 
sity  functions  are  quite  different.  The  poisson  model  is  often  used  to  represent 
target  evasive  or  "jinking”  maneuvers.  The  poisson  square  wave  can  only  take 
on  values  of  ±  so  at  any  given  time  its  probability  density  function  consists 
of  impulses  with  a  weighting  of  0. 5  at  plus  and  minus  cr^,  whereas  the  above 
markov  process  is  assumed  to  have  a  gaussian  amplitude  distribution.  There¬ 
fore,  the  response  of  an  amplitude  dependent  nonlinear  operator  could  be  quite 
different  when  driven  by  each  of  these  two  signal  forms.  However,  if  the  ran¬ 
dom  square  wave  is  passed  through  a  narrow  band  filter  or  integrator,  it  would 
experience  broadening  due  to  the  filter's  finite  bandwidth.  In  the  case  of  an 
integrator,  the  resulting  wave  shape  would  be  a  series  of  linear  segments  of 
constant -Slope.  By  application  of  the  central  limit  theorem,  as  discussed  in 
Ref.  9,  the  output  distribution  approaches  the  gaussian  density  function  as  the 
number  of  stages  of  filtering  increases.  Thus,  although  the  poisson  square  wave 
may  in  some  respects  be  a  more  realistic  target  maneuver  model,  we  take 
advantage  of  the  statistical  similarity  of  these  processes  to  justify  repre¬ 
senting  this  random  effect  by  a  band-limited  gaussian  process  to  facilitate  the 
subsequent  CADET  analysis. 

A. 4  THE  AUTOPILOT -AIRFRAME  MODEL 

In  accordance  with  the  assumption  that  the  missile  and  target  trajec¬ 
tories  are  confined  to  a  plane,  we  describe  the  missile  airframe  orientation  by 
the  variables  depicted  in  Figure  A.  4-1.  This  figure  establishes  the  sign  con¬ 
vention  of  each  quantity;  each  variable  is  positive  as  shown.  Note  that  we  are 
particularizing  the  airframe  model  at  this  point  by  discussing  the  tail-controlled 
tactical  missile.  This  is  done  to  provide  a  concrete  model  for  study,  and  not 
necessarily  to  exclude  other  configurations. 
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y-  AXIS 
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A.  4-1  Geometric  Definition  of  Intercept 
Plane  System  Variables 


The  airframe  variables  under  consideration  are 
x .  =  0  =  body  angular  rate 

7  m 

Xg  =  a  =  angle  of  attack 

Xg  =  6  =  control  surface  deflection 

a  =  missile  lateral  acceleration  (normal  to 
m  the  missile  velocity  vector) 

We  neglect  gravity  effects,  tacitly  assuming  that  the  intercept  plane 
is  horizontal  or  that  the  missile  has  perfect  gravity  compensation.  In  a  gen¬ 
eral  situation,  the  differential  equations  expressing  the  airframe  dynamics 
are  nonlinear  and  time -varying  due  to  the  dependence  of  the  airframe  para¬ 
meters  on  altitude,  angle  of  attack,  and  mach  number.  However,  in  this 
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study  we  express  the  state  variable  differential,  equations  in  linear  form*, 


x.  =  M  x.  +  M  x_  +  M.x* 
4  q  4  cl  5  5  6 


xc  =  x.  -  L  xc  -  L,xe 
5  4  a  5  5  6 


xg  =  -  UXg  +  n u (t) 


where  u(t)  represents  a  commanded  fin  deflection,  1/m  is  the  control  fin  actua¬ 
tor  time -constant,  and  the  constants  Ma,  Mg,  and  Lg  represent  the 
airframe  stability  derivatives.  The  latter  are  obtained  from  the  nonlinear  air¬ 
frame  parameters  by  making  the  following  assumptions: 

•  Missile  velocity  is  constant  (drag  effects  are 
negligible  over  the  period  of  time  considered). 

•  Altitude  remains  nearly  constant. 

•  The  center  of  pressure,  mass  and  inertia  of  the 
missile  are  constant. 

•  Lift  force  is  linearly  related  to  changes  in  angle 
of  attack  about  some  trim  condition  and  to  con¬ 
trol  fin  deflection. 

•  Control  fin  actuator  dynamics  are  linear. 

•  Fin  effectiveness  is  independent  of  angle  of  attack. 

Practical  experience  has  shown  that  the  resulting  autopilot  response  character¬ 
istics  closely  approximate  those  of  the  nonlinear  airframe  near  the  given 
nominal  conditions.  The  output  of  the  airframe  modei  is  the  missile  lateral 
acceleration,  which  is  given  by 


A  highly  nonlinear  airframe  including  drag  effects  has  been  studied  via 
CADET  in  Ref.  3. 
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a  =  v  [x.  -  xJ  =  v  [Lx-  +  L»xft]  (A.  4-2) 

m  mL  4  5J  mL  a  5  o  6J 

where  vm  is  the  magnitude  of  the  missile  velocity  vector. 

For  typical  values  of  the  stability  derivatives  in  Eq.  (A.  4-1),  the 
missile  airframe  will  exhibit  an  underdamped  or  even  an  unstable  response  to 
a  commanded  fin  deflection.  Acceptable  control  is  achieved  by  introducing 
feedback  compensation  in  the  fin  deflection  command, 

u(t)  =  -  [k  a  -  k  (a  /v  )  -  k.  6  ] 
w  L  c  c  a  m  nr  d  mJ 

(A.  4-3) 

=  -fka  -  k  (L  xc  +  L.xJ  -  k.  x  J 
Lcc  a  a  5  5  6'  d4j 

where  a  is  the  commanded  acceleration  provided  by  the  guidance  law  model 
(Section  A.  5).  The  parameter  k  is  chosen  to  give  unity  steady  state  gain 
from  ac  to  am,  and  k^  and  ka  are  chosen  to  give  the  desired  transient  response 
A  complete  block  diagram  of  the  compensated  missile  dynamic  equations  is 
shown  in  Fig.  A.  4-2. 


M1SV7 


Figure  A.  4-2  Compensated  Missile  Airframe  Dynamics 
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For  ready  assessment  of  the  compensated  missile  airframe  dynamics, 
it  is  convenient  to  use  a  transfer  function  formulation  of  the  model.  Given  two 
outputs,  am  and  flm,  we  desire  to  obtain  g^(s)  and  g2(s)  to  provide  the  input- 
output  relations  indicated  in  Fig.  A.  4-3.  Utilization  of  standard  block  diagram 


R-11598 


Figure  A.  4-3  Transfer  Function  Formulation  of  the 

Compensated  Missile  Airframe  Dynamics 


reduction  techniques  shows  that  the  dynamics  indicated  in  Fig.  A.  4-2  are  equiv¬ 
alent  to  the  transfer  function  formulation  depicted  in  Fig.  A.  4-3,  where 


with 


g,(s)  = 

g2(s)  = 


e3s  +  e2s  +  el 

"3  2 

S  +  CgS  +  CgS  +  c 

d2s  +dl 
~3  2 

s  +  CgS  +  CgS  +  C 


1 


1 


c1=m[v  +  <ka  +  1^)  B] 

c2  =  y  +  (.[t  +  L6Mqka-M6kb] 

°3  =  C  +  Ml  -  I-g  ka) 


(A. 4-5) 

(A.  4-6) 


(A.  4-7) 


A -10 


THE  ANALYTIC  SCIENCES  CORPORATION 


dl  =  <**kc 

(A'4'8) 

ex  =  **vm 

e2  “  “vmL6M<!kc 

e„  =-  fiv  L.k  (A.  4-9) 

3  m  o  c 

P  =  L.M  -  L  M. 
o  a  a  o 

y  =  “Ma  -  LaMq  (A.  4-10) 


The  aerodynamic  coefficients  used  in  this  study  are  chosen  to  corres¬ 
pond  to  an  intercept  at  35,000  ft. ,  with  a  missile  velocity  magnitude  vm=3000 
ft/sec,  the  data  being  taken  from  Ref.  14,  Vol.  n,  Appendix  H.  The  compensa¬ 
ting  gains  ka,  k^  and  kfi  Eq.  (A.  4-3)  are  set  to  achieve  a  dynamic  performance 
that  is  typical  of  a  missile  in  the  terminal  homing  phase.  These  parameters  and 
the  corresponding  transfer  function  coefficients  are  given  in  Table  A.  4-1. 


The  dominant  time  constant  of  the  compensated  airframe  is  Tj  =  l/s^= 
0. 316  sec.  The  fact  that  e^,  eg  and  e^  do  not  all  have  the  same  algebraic  sign 
demonstrates  that  g^(s)  has  a  right  half  plane  zero,  which  is  characteristic  of 
the  tail-controlled  missile  configuration  depicted  in  Fig.  A.  4-1. 


TABLE  A. 4-1 


COMPENSATED  MISSILE  AIRFRAME  DATA  IN 
THE  TERMINAL  HOMING  PHASE 


Parameter 


Symbol 


Value 


Actuator  Lag 
Time  Constant 


Aerodynamic 

Coefficients 


Compensating 

Gains 


Transfer 

Function 

Coefficients 


0.0533  sec 


-0.462  sec"1 
-5.81  sec “2 
-72.0  sec"2 
0. 379  sec"* 
0.070  see-* 


1.02  sec 
0.188  sec 

0.476  xlO-2  sec2/ft 


720.0  sec"3 
275.  3  sec"2 
18.  3  sec"* 

0.240  sec-2  ft"* 
0.642  sec"1  ft"1 
720.0  sec"2 
-0.865  sec"2 
-1.87  sec"1 


-3.16  sec"1 
-7. 56  +13. Oj  sec"1 
-7.  56-13.0j  sec"1 


Transfer 

Function 

Poles 
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A.  5  THE  GUIDANCE  SUBSYSTEM  MODEL 

The  guidance  signal  from  the  seeker  (discussed  in  the  following  section) 
is  passed  through  a  single-stage  low-pass  noise  filter,  the  output  of  which  is 
thus  a  filtered  estimate  of  LOS  angle  rate,  x^  ^  0.  The  classical  proportional 
guidance  law  is  then  implemented,  which  calls  for  a  commanded  acceleration 

whose  component  normal  to  the  line-of-sight  (LOS)  is  proportional  to  the  clos- 

.4 

ing  velocity  times  the  estimated  LOS  angle  rate,  =  9;  that  is,  we  desire  to 

cause  a_  to  satisfy 
m 

am  cos  (0X  "  0)  =  n\B  (A.  5-1) 

where  the  parameter  n'  is  designated  the  navigation  ratio.  The  closing  velocity 
is  obtained  by  projecting  the  missile  and  target  velocity  vectors  onto  the  instan¬ 
taneous  line  of  sight;  as  shown  in  Fig.  A.  2-1, 

vc  =  vm  cos  (^  -  0)  +  vtcos  (ea  +  0)  (A.  5-2) 


In  order  to  achieve  a  response  that  obeys  Eq.  (A.  5-1),  the  ideal  acceleration 
command  a '  should  be  chosen  to  satisfy 

v 


(A. 5-3) 


where  the  incorporation  of  the  factor  1/cos  (0^- 0)is  known  as  secant  com- 
pensation. 


In  mechanizing  the  guidance  law,  the  value  of  the  closing  velocity  is 
never  known  exactly.  If  a  radar  homing  seeker  is  used,  then  a  reasonable 
estimate  of  v  can  be  obtained  by  doppler  measurements  or  by  differencing 

V 

range  measurements.  An  infrared  seeker  system  generally  does  not  yield  a 
good  estimate  of  range,  in  which  case  v  may  be  taken  to  be  a  prespecified 

V 
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constant.  Any  uncertainty  in  the  closing  velocity  is  modelled  by  introducing 
a  variable  ey  into  Eq.  (A.  5-3)  which  represents  either  a  band  limited 
noise,  obtained  by  a  single-stage  low-pass  filter  with  white  noise  input,  or 
a  bias,  denoted  simply  e^.  Thus 


'  =  n '  0  [v  +  v 
c  [mi 


COSfrg  +  0) 
t  cos(x8  -  0) 


+  e 


] 


(A.  5-4) 


provides  the  final  ideal  acceleration  command  used  in  this  study,  where 

A 

Xjg  =  ©v  is  modelled  by  one  of  the  differential  equations 

Random  Uncertainty:  x^  =  -  03 +  w4>  Efx^O)]  =  ® 

(A.  5-5) 


Bias  Uncertainty: 


*12  =  °>  x12(0)  =  evb 


and  w4  is  a  white  noise.  Thus  we  can  study  either  the  effect  of  the  noisy  esti¬ 
mation  of  v  or  of  a  constant  error  in  the  assumed  value  of  v  .  i 
c  c 


Finally,  the  guidance  law  must  account  for  an  important  ponlinear  con¬ 
straint  on  missile  operation  --  namely,  acceleration  command  limiting.  The 
actual  acceleration  command  a  that  determines  the  input  to  the  fin  deflection 
actuator  in  Fig.  A.  4-2  must  not  exceed  the  structural  capacity  of  the  airframe 
and  must  not  be  so  large  as  to  cause  the  missile  to  stall.  Thus  the  ideal  accel¬ 
eration  command  a^  must  be  limited  in  order  to  prevent  excessive  lateral 
acceleration  command  levels;  the  limiting  procedure  is  represented  by  the  sat¬ 
uration  nonlinearity 


a  = 
c 


<  a 

c  '  !  c  I  “  max 

(A.  5-6) 

W  sign  (ac>>  Kl>amas 
The  guidance  law  features  described  above  are  all  incorporated  in  the  system 
model,  as  illustrated  by  the  block  diagram  in  Fig.  A.  5-1. 
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CLOSING  VELOCITY 
ERROR  MODEL 

, - * - i 


R-J1599o. 


NOISY 
LOS  RATE 
ESTIMATE 


Y 

NOISE  FILTER 


CLOSING  VELOCITY  WITH  SECANT  ACCELERATION 

COMPENSATION  AND  UNCERTAINTY  COMMANDS 


Figure  A.  5-1  Guidance  Law  Model 


A.  6  THE  SEEKER  SUBSYSTEM  MODEL 


There  are  several  effects  inherent  to  the  seeker  which  can  have  a 
marked  influence  on  overall  missile  performance.  These  include 


•  Boresight  error  distortion 

Noise 

Aberration 

Receiver  and  signal  processing  characteristics 

•  Disturbance  torque  inputs 

Seeker  mass  imbalance 

Seeker  gimbal  friction 

Spring  restoring  forces  on  the  seeker  head 

Boresight  Error  Distortion  -  A  fundamental  variable  in  the  seeker  sub¬ 
system  is  the  true  boresight  error,  ^^rue  defined  by  the  angle  between  the  anten¬ 
nal  centerline  and  the  instantaneous  line-of -sight  (LOS);  referring  to  Fig.  A.  6-1 


y-AXIS 


*-11600 


MISSILE 

CENTER  OF  GRAVITY 

Figure  A.  6-1  Seeker  System  Configuration 

The  estimated  or  measured  value  of  the  boresight  error  will  differ  from  *true 
due  to  several  factors;  among  the  more  important  of  these  are  aberration,  noise, 
and  nonlinear  receiver  characteristics. 


The  effect  of  aberration  is  very  highly  dependent  upon  the  geometry  of 
the  seeker -detector  cover,  the  frequency  and  polarization  of  the  incident  energy 
and  other  factors;  furthermore,  it  is  variable  due  to  manufacturing  tolerances, 
possible  erosion  during  flight  and  changes  in  environmental  parameters.  This 
phenomenon  can  be  represented  by  a  nonlinear  and  possibly  time-varying  opera¬ 
tion  on  the  look  angle,  ®look  =  ®  so  that  an  effective  boresight  error, 

is  obtained  in  the  form 

ceff  =  0look  +  0ab  "  \ 

where 

®ab  -  ^ab^look^ 

as  depicted  in  Fig.  A.  6-2. 
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R-J1601 


LOS 

ANGIE' 


.  ^look 

1  W 

ABERRATION 

6look+eob  , 

)  ► 

CHARACTERISTIC 

Kv 

MISSILE 
BODY  ANGLE 


f‘ 


EFFECTIVE 

BORESIGHT 

ERROR 


SEEKER 

ANGLE 


Figure  A.  6-2  Boresight  Aberration  Model 


The  effect  of  a  linear  aberration  characteristic  (0^  =  kr  9look)  was  studied  in 
Ref.  2.  The  case  in  which  a  radar  homing  system  with  nonlinear  aberration  is 
considered  (where  the  aberration  of  incident  energy  is  caused  by  the  protective 
radome)  was  successfully  analyzed  via  CADET  in  Ref.  3.  In  the  latter  study,  the 
radome  aberration  characteristic  was  modelled  as  a  piece -wise  linear  relation 
with  odd  symmetry  and  five  distinct  linear  segments,  as  depicted  in  Fig.  A.  6-3. 
To  avoid  duplication,  this  effect  is  not  considered  further  in  this  study. 


R-9171 


Figure  A.  6-3  Nonlinear  Angular  Aberration  Characteristic 
Investigated  in  Ref.  3 
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In  considering  the  degrading  effects  of  noise,  we  include  three  funda¬ 
mental  categories  of  effects.  Inverse  range  proportional  noise,  which  has  an 
effective  rms  level  of  the  form 


is  representative  of  a  noise  source  that  increases  in  effect  as  range  approaches 
zero.  Such  a  phenomenon  is  target  angular  scintillation  (caused  by  the  apparent 
motion  of  the  target  due  to  the  change  in  position  of  the  target  centroid  of  radia¬ 
tion).  This  can  be  modelled  as  a  wide -band  noise  state,  x^,  with  constant  rms 
level,  cTj^,  multiplied  by  a  gain  1/r.  Range  proportional  noise  includes  any 
noise  source  that  remains  constant  in  rms  level  in  such  a  way  that  the  effective 
signal-to-noise  ratio  decreases  as  the  missile  approaches  the  target,  i.  e. ,  as 
range  approaches  zero.  This  type  of  random  disturbance  is  represented  by  an 
equivalent  noise  with  an  rms  level  of  the  form 

ab  =  ra15 

which  in  turn  can  be  modelled  by  a  wide-band  noise  state  with  a  constant 
rms  level  of  passing  through  a  gain,  r.  Noise  sources  that  exhibit  this  pro¬ 
perty  are  the  distant  stand-off  jammer  and  receiver  noise  (generally  due  to 
thermal  effects).  Range  independent  noise  represents  noise  sources  that  have 
a  constant  effect  on  the  signal-to-noise  ratio;  target  amplitude  scintillation  (due 
to  time -varying  effective  target  cross  section,  for  example)  and  seeker  servo 

noise  are  typical  examples  of  such  noise  sources  that  can  be  modelled  by  a 

2 

noise  state  x^g  of  constant  variance  cr^.  The  complete  noise  model  is  shown 
in  Fig.  A.  6-4  where  w^,  Wg,  and  w^  are  gaussian  white  noise  processes. 

All  three  types  of  noise  considered  above  have  been  treated  in  previous 
studies  (Refs.  2,  3).  It  should  be  noted,  however,  that  the  earlier  implementa¬ 
tions  of  this  model  were  linear  time -varying;  i.  e. ,  r(t)  was  assumed  to  be  known 


A-I8 


Figure  A.  6-4  General  Seeker  Noise  Model 

exactly  in  the  noise  model.  Here  we  rigorously  implement  the  nonlinear  rela¬ 
tion  indicated  in  Fig.  A.  6-4, 


ns  =  x16  +  x15  ^X102  +  X11  +  X14'  \o  +  xll2  (A*  6'2) 

The  receiver  characteristic  is  a  potentially  complicated  effect,  highly 
dependent  upon  the  specific  antenna  design,  type  of  detector,  and  signal  pro¬ 
cessing  scheme.  In  order  to  avoid  a  very  specialized  model  based  on  a  partic¬ 
ular  tactical  missile,  we  confine  our  attention  to  one  phenomenon:  the  attenuation 
of  the  received  signal  which  occurs  ”’hen  boresight  error  €  becomes  large,  i.  e. , 
when  c  approaches  e  in  Fig.  A.  6-5a.  The  detector  alone  will  have  an  output 
which  is  very  nearly  proportional  to  its  input  for  small  values  of  e;  however,  as 
the  boresight  error  approaches  c  ,  we  note  in  Fig.  A.  6 -5b  that  the  signal 
strength  decreases  to  a  null.  If  the  antenna  pattern  has  appreciable  sidelobe  sen¬ 
sitivity,  there  may  also  be  some  response  for  values  of  c  greater  than  fmax- 
The  upper  limit  on  the  boresight  error,  7,  such  that  the  detector  characteristic 
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(b)  ANTENNA- DETECTOR 
CHARACTERISTIC 


(c)  SIGNAL  PROCESSING 
NONLINEARITY 


Figure  A.  6-5  Receiver  Boresight  Err^r  Distortion  Effects 


is  nearly  linear  for  |<  |  less  than  c  is  quite  variable,  depending  on  the  type  of 
target  tracking  system  under  consideration.  For  monopulse  radar  or  infrared 
detectors,  c  could  be  as  small  as  a  fraction  of  a  degree. 


The  undesirable  detector  null  and  possible  sidelobe  response  can  be 
circumvented  in  the  signal  processing  scheme.  In  the  caeo  considered  here, 
some  value  €^m<€  is  chosen;  a  nonlinearity  is  then  introduced  such  that  whenever 
the  boresight  error  magnitude  exceeds  e ^m,  the  output  of  the  signal  processor 
is  held  at  ±c^m*  This  provides  a  simple  model,  depicted  in  Fig.  A.  6 -5c,  which 
will  capture  the  effect  of  a  narrow  antenna  beamwidth  and  a  reasonable  signal 
processing  nonlinearity. 
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The  combined  effects  of  noise  and  receiver/signal  processing  charac¬ 
teristics  are  illustrated  in  the  general  boresight  error  monel  (without  aberration) 
shown  in  Fig.  A.  6-6.  We  mention  in  passing  that  a  more  exact  noise  model 
might  divide  noise  sources  into  external  and  internal  effects,  i.  e. ,  noise  sources 
entering  the  boresight  error  model  just  before,  as  well  as  after,  the  receiver 
characteristic.  However,  this  categorization  is  somewhat  artifically  detailed 
when  compared  with  a  realistic  situation  where  noise  levels  are  only  approxi¬ 
mately  known.  While  the  aberration  and  signal  processing  nonlinearities  may 
have  a  significant  impact  on  the  closed  loop  dynamics  of  the  overall  missile  - 
target  intercept  model,  the  actual  noise  input  point  is  not  particularly  critical. 

In  fact,  injecting  all  noise  sources  after  the  receiver  in  Fig.  A.  6-6  results  in 
a  worst -case  model  of  the  effect  of  the  signal  processing  nonlinearity,  since  it 
neglects  the  attenuation  of  external  noise  sources. 


R-11605 


Figure  A.  6-6  Final  Boresight  Error  Measurement  Model 
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Disturbance  and  Control  Torques  -  The  seeker  model  is  completed  by 
developing  a  suitable  tracking  and  stabilization  control  system  including  several 
important  sources  of  disturbance  torque  inputs.  In  terms  of  the  inertially- 
referenced  angles  <t>  and  ip  shown  in  Fig.  A.  6-1,  we  will  derive  a  relation  of  the 
form 

V^V'VTd  (A.  6-3) 


where  1^  is  the  moment  of  inertia  of  the  seeker  head  about  the  gimbal  pivot,  Tc 
is  the  external  control  torque  (derived  from  an  electric  servo  motor,  for  exam¬ 
ple),  and  Tj  is  the  total  disturbance  torque.  The  latter  has  three  components, 


+  Tf  +T 


r 


(A.  6-4) 


where  Tm  is  an  effective  torque  due  to  mass  imbalance,  and  we  consider  two 
external  torque  components,  T^  due  to  nonlinear  friction  in  the  gimbal  and  Tf 
due  to  nonlinear  restoring  torques.  Since  the  seeker  head  center  of  gravity 
(eg)  is  generally  displaced  from  the  pivot  point,  as  shown  in  Fig.  A.  6-1  and 
specified  by  the  parameters  r^  and  0q,  the  moment  of  inertia  I  is  related  to 
the  corresponding  moment  inertia  referred  to  the  eg  by 

Ip=I0  +  mro  (A.  6-5) 


where  m  is  the  mass  of  the  seeker  head. 

The  external  torques  due  to  spring  and  friction  effects  are  modelled  by 
the  relations 


Restoring  Torque:  Tr  =  1^(0^) 
Friction  Torque:  Tf  =  f2(0h) 


(A.  6 -6) 
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where  0^  is  the  angle  between  the  seeker  and  missile  center  lines.  Often 
restoring  torques  are  linear  for  small  angle  deflections,  becoming  nonlinear 
only  as  0^  increases  in  magnitude,  as  illustrated  in  Fig.  A.  6 -7a.  This  be¬ 
havior  corresponds  to  the  symmetric  "hard  spring”  case  (Ref.  15)  where  the 
elastic  limit  of  a  spring  is  exceeded  and  Hooke's  law  for  linear  spring  behavior 
becomes  invalid;  often  the  nonlinear  term  is  taken  to  be  a  power  law  relation, 


W=ki 


lim 


K 

sign  (6h) 


where  k  is  an  integer  greater  than  one,  so  that  exhibits  a  distinct  departure 
from  linearity  as  1 0^  I  exceeds  0^.  m  which  is  typical  of  a  symmetric  nonlinear 
spring  characteristic.  A  common  type  of  nonlinear  friction  is  the  dry  or  Cou¬ 
lomb  effect  (Ref.  15),  where 

fA>  =  k2  slgn 

i.  e. ,  the  friction  term  of  the  disturbance  torque  has  constant  magnitude  with  the 
algebraic  sign  of  the  gimbal  angle  rate.  Illustrations  of  these  disturbance  torque 
terms  are  depicted  in  Fig.  A.  6-7. 


R-11604 


(a)  Nonlinear  Restoring  Torque  (b)  Nonlinoor  Friction  Effocl 

Figure  A.  6  -7  External  Disturbance  Torque  Models 
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The  effective  disturbance  torque  component  due  to  seeker  mass  im¬ 
balance  can  readily  be  determined  by  application  of  the  basic  principles  of 
mechanics  (Ref.  13).  The  equations  of  motion  of  the  seeker  head  can  be  de¬ 
rived  directly  by  applying  the  Euler  equation 


d_ 

dt 


/Vu\  3U 
\a$  /  W 


T  -  T  -  T, 
erf 


(A.  6-7) 


where  U  is  the  seeker  head  kinetic  energy,  and  the  right  side  represents  all 
external  torque  sources.  The  kinetic  energy  of  the  seeker  head  has  two  terms: 


TT  1  2  1 T  ,h2 
U  =  2mv  +2I0(W 


(A.  6-8) 


where  v  is  the  magnitude  of  the  velocity  of  the  seeker  head  center  of  gravity 
and  0  is  the  angular  velocity  of  the  seeker  head  rotation  about  the  eg.  The 
velocity  magnitude  is  obtained  by  resolving  the  missile  velocity  vector  v  m 
(describing  the  motion  of  the  missile  eg)  into  x  and  y  components,  then  adding 
to  this  the  relative  motion  of  the  seeker  eg  with  respect  to  the  missile  eg. 
Referring  to  Fig.  A.  6-1,  the  seeker  head  eg  has  the  relative  coordinates 


x  =  r1cos  (0m  +  9  )  +  rQcos  0 


y  =  r.  sin  (6  +  0  )  +  rn  sin  0 

J  1  '  m  p'  0 


Generally  0^  (the  fixed  angle  between  the  missile  centerline  and  the  seeker 

gimbal  pivot)  is  quite  small,  so  it  will  be  neglected  here.  The  total  velocity- 

components  of  the  seeker  head  eg  are  then 

dx  •  *  *  • 

vx  =  vm  cos  9 «  +  dF  =  Tm  008  9*  '  rl  9m  SIn  9m  '  r0 I*  Sln  * 

VY  =  vm  sin  h  +  §=  vm  sin  6i  +  rl  9m  cos  9m  +  r0  *  cos  * 
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Using  standard  trigonometric  identities,  the  squared  velocity  magnitude  is 
then  determined  to  he 


+  <rl  %)2  +  {l0i)2  +  2  r0  rl  V cos  (l>  •  em> 


+  2rOvm*sin(0« '*>  +  2rlvmVin(W  (A'6-9) 

The  equation  of  motion  can  then  be  derived  directly  from  Eqs.  (A.  6-7)  to 
(A.  6-9);  combining  these  relations,  we  obtain  the  final  differential  equation 
of  motion, 


-  'p  *  = m  Vi 


[0  cos 
L  m 


(* 


0  ) 
m 


+ 


(0m)2sin  (if)' 


e  )] 

m 


+  mrAv  0.  cos  (0-0)  (A.  6-10) 

0  m  £  v  l  ' 

Comparing  Eq.  (A.  6-10)  with  Eqs.  (A.  6-3)  and  (A.  6-4),  we  identify  the  seeker 
mass  imbalance  disturbance  torque,  Tm  in  Eq.  (A.  6-4),  to  be  given  by  the 
terms  on  the  right  side  of  Eq.  (A.  6-10).  Combining  the  latter  with  the  friction 
and  spring  disturbance  torque  components,  we  obtain 

T,  =  f1  (0. )  +  f„  (0, )  +  mr.  [ri  0  cos (0-0  ) 
d  1  v  h'  2  '  h7  0 L  1  m  v  nr 


+  ri  (0  )2  sin  (0-0  )  +  v  0,  cos  (0  -  0, )] 
1 '  m'  v  m'  m  jt  v  Z /J 


(A.  6-11) 


The  control  torque  T  in  Eq.  (A.  6-3)  is  chosen  to  make  the  seeker 

c 

track  the  target,  i.  e. ,  to  maintain  the  measured  boresight  error  at  a  small 
value.  The  nominal  seeker  is  designed  under  the  assumption  that  there  is  no 
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friction  and  that  spring  effects  are  negligible;  thus,  it  is  necessary  to  include 
rate  feedback  in  the  torque  command  (a  feedback  term  proportional  to  0  which 
is  measured  by  a  rate  gyro)  to  provide  suitable  damping.  Thus  we  write  the 
nominal  control  torque  as 


where  is  the  track  loop  time  constant,  is  the  rate  gyro  gain,  and  kg  is 
the  torque  servo  gain.  The  implementation  of  this  control  law  is  depicted  in 
Fig.  A.  6-8. 


R-11606 


Figure  A.  6-8  Nominal  Seeker  Track  Loop 

(Neglecting  All  Nonlinear  Effects) 

While  the  implementation  of  the  seeker  control  function  depicted  in  Fig. 
A.  6-8  will  provide  an  adequate  response  under  ideal  conditions,  it  can  be  shown 
(see  Section  A.  7)  that  the  dynamic  response  of  the  seeker  is  quite  sensitive  to 
steady  state  disturbance  torque  inputs.  Since  we  have  already  indicated  that  dis¬ 
turbance  torques  are  to  be  considered  in  our  study,  compensation  of  the  seeker 
must  be  included  to  achieve  satisfactory  performance.  A  simple  and  effective 
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compensation  procedure  is  to  insert  proportional  plus  integral  cascade  compen¬ 
sation  before  the  torque  summing  junction  in  Fig.  A.  6-8.  That  is,  we  specify  the 
compensated  control  torque  by 


s  +  k 


0 


(A. 6-13) 


This  relation  corresponds  to  the  differential  equation 


=  k.T 
0  cn 


+  T 


cn 


(A.  6-14) 


The  complete  seeker  simulation  model,  representing  the  synthesis  of 
the  dynamic  equations  derived  in  this  section,  is  shown  in  Fig.  A.  6-9.  We 
demonstrate  in  Section  A.  7  that  the  indicated  output,  n,  is  a  noisy  estimate  of 
LOS  rate  (0). 


A.  7  TRANSFER  FUNCTION  REPRESENTATION  OF  THE  EQUIVALENT 
LINEAR  SEEKER 

For  a  subsystem  of  the  complexity  of  the  seeker  as  modelled  in  Fig. 
A.  6-9,  it  is  often  helpful  to  derive  the  transfer  function  formulation  of  the 
linear  system  obtained  by  neglecting  all  nonlinearities.  Several  assertions 
made  in  simplifying  the  seeker  model  are  based  on  this  representation,  and 
the  procedure  used  for  the  purpose  of  designing  the  compensation  network 
(choice  of  k^)  can  best  be  treated  in  this  way. 

We  define  four  inputs  (refer  to  Fig.  A.  7-1), 

0 

•  « 

0 

m 

n  (t) 
s'  ' 

T  (A.  7-1) 

m 
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The  transfer  functions  d^(s)  to  d^(s)  can  be  found  for  the  equivalent  block  dia¬ 
gram  representation  depicted  in  Fig.  A.  7-2  using  standard  state  vector  formu¬ 
lation  techniques, 


where 


x  =  Fx  +  Gu  \ 

D(s)  =  E  +H(sI-F)~1G  (A.  7-2) 

tj  =  Hx  +  Eu  I 


U1 

u2 

Uo 


D(s)  =  [d1(s),  d2(s),  dg(s),  d4(s)]  (A.  7-3) 


and  d.(s),  the  elements  of  D,  correspond  to  the  scalar  transfer  function  indica¬ 
ted  in  Fig.  A.  7-2.  The  final  results  of  this  analysis  are  summarized  in  Eqs. 
(A.  7-4)  and  (A.  7-5). 

R-ll  610 


Figure  A.  7  -2  Linear  Seeker  Model  in 
Transfer  Function  Form 
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J(„,  k3  s[s  +q3s+Pl] 
d  1  (s)  =  -=~  -j - 5 - 

1  T1  s%  q/  +  q2s  +  qj 


-k. 


k2s  +  ki 


d2(S)  =  ^ - 5--~  o  1 - 

p  1  s[s 6  +  q3s2  +  q2s  +  q^ 


d3(S>  =  dl(S> 


‘Vs)  I  Ti  ”3  3 

p  1  s  +q3s  +  <j2s  +  qi 


k.  +k  k  kn 

=  j _ s  g  ° 

i 


(A.  7-4) 


«1 


ks  k0  k3 
JpTl 


ksk3 

q2  =  P1  +  "T~t7 
P  1 


(A.  7-5) 


kn  +  k  k 

q  =J _ 

q3  I 

P 


The  nominal  seeker  is  defined  by  a  choice  of  parameters  that  leads  to 
acceptable  dynamic  behavior  in  the  absence  of  disturbance  torques;  the  data  and 
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transfer  functions  are  given  in  Table  A.  7  -1.  A  second  case  is  the  nominal 
compensated  seeker,  -which  has  been  designed  to  exhibit  a  significantly 
better  performance  in  the  presence  of  disturbance  torques;  the  design  by 
root  locus  techniques  is  indicated  subsequently  and  the  transfer  functions  are 
summarized  in  Table  A.  7  -2. 

hi  both  the  nominal  and  the  nominal  compensated  seeker,  we  note 
that  d2(s)  H  0.  This  demonstrates  that  with  no  linear  friction  or  spring  re¬ 
storing  torques,  the  seeker  has  perfect  stabilization,  i.  e. ,  the  measured 
boresight  error  is  unaffected  by  rotation  of  the  missile  body. 

For  frequencies  considerably  less  than  10  rad/sec,  we  have  d^  =  dc^=o, 
which  is  the  transfer  function  of  a  differentiator.  Hence,  the  assertion  that  jp 
is  an  estimate  of  the  LOS  angular  rate  (8)  holds  at  low  frequencies. 

The  seeker  compensation  removes  steady  state  disturbance  torque  sen¬ 
sitivity,  as  is  shown  by  the  zero  of  dc4(s)  at  s  =  0.  To  demonstrate  this  point  in 

more  detail,  we  assume  a  step  disturbance  torque  input, 

u4(t)=T0,  t*0,  (A.  7-6) 

yielding  a  response  which  is  represented  in  Laplace  transform  notation  by 

T 

V(s)  =  d4(s) 


or,  by  substitution, 


V  (s) 
-n'  ' 


1000  x0 
12  s  (s  +  10)(s  +  50) 


Vs)  =  ~°r 


(A.  7-7) 


0 


s3  +  60  s2  +  1700  s  +  10,000 
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k3  =  1 


Transfer  Functions 


j  f  \  _  100  s  (s  +  60) 
nV'  '  ~  12  (s  +  10)(s  +  ! 


k  =  1 
g 


dn2(S)  -  0 


J  d  ,(s)  -  d„,(s) 


I  =  0.1  in-oz-sec^ 
P  _ 

T.  =  0. 12  sec 


*  _  1000  1 
n4^s'  12  (s  +  10)(s  +  5 


TABLE  A.  7-2 

THE  NOMINAL  COMPENSATED  SEEKER 


Parameters 


Transfer  Functions 


kl  *  k2  =  0 


k  =  k  =1 

3  g 


J  /  \  100  s  s“  +  60s  +  1200 

dci(s)  =  -To  “S — — o  - - 

C1  1  s  +  60s^  +  1700s  +  10, 000 


k^  =  20  sec 


I  =0.1  in-oz-sec^l 
P 


=  0. 12  sec 


dc2(s)  =  0 


dc3(s)  =  dcl(s) 


,  /v  1000  s 

a«4's'  “  12  3  2 

s°  +  60s  +  1700  s  +  10, 000 


POLES:  s  =  -7. 71,  s  =  -26.  1  -  24.  8j 
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To  obtain  the  effect  of  the  disturbance  torque  after  transients  disappear,  we 
evaluate  the  steady  state  value  of  rj  by  the  final  value  theorem  of  Laplace  trans¬ 
form  theory:  if  f(s)  is  the  Laplace  transform  of  f(t)  and  if  a  finite  steady  state 
value  of  f(t)  exists,  then 

fss  "  Hm  fft'  =  f s  f(s)] 
t-»  »  s  -+0 


The  application  of  this  result  to  Eq.  (A.  7-7)  yields 

_Te 

^nss  6 


n 


css 


=  0 


(A.  7-8) 


which  establishes  that  the  proportional  plus  integral  compensation  technique 
(insertion  of  the  transfer  function 


h  (S)  =  t±lo 

cv  s 

in  Eq.  (A.  6-13))  is  effective  in  eliminating  the  response  to  constant  disturbance 
torque  inputs. 

To  complete  the  seeker  compensation  design,  we  note  that  for  the  nom¬ 
inal  values  of  the  seeker  variables  given  in  Table  A.  7-1  (except  ^  0),  the 
denominator  of  the  transfer  functions  d^(s)  is  given  by 

3  2  3  2 

s  +  qgS  +  qgS  +  q^  =  s  +  60s  +  (500  +  60  k^)  s  +  500  k^ 

=  (s3  +  60s2  +  500s)  +  60  kQ  (s  +  8.  333) 
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In  choosing  k^,  it  is  desired  to  achieve  a  dominant  pole  at  s  =  -  10  and  to  have 
the  remaining  subsidiary  poles  as  far  from  the  origin  of  the  s -plane  as  possible 
subject  to  the  condition  that  they  must  be  suitably  damped.  A  root  locus  analy¬ 
sis  is  effective  in  finding  a  suitable  value  of  satisfying  these  conditions.  In 
Fig.  A.  7-3,  we  see  that  as  k^  *♦ 00 ,  the  dominant  pole  approaches  -8.  33;  how¬ 
ever,  for  kQ  >  20,  we  have  subsidiary  poles  with  unsatisfactory  damping,  i.  e. , 
C  <  0. 707.  The  condition  C  =  0. 707  is  specified  in  order  to  ensure  a  transient 
step  response  with  no  overshoot;  thus  our  final  compensated  seeker  design  is 
achieved  by  choosing  kQ  =  20,  as  indicated  in  Table  A.  7 -2. 

*-11607 

IMAGINARY 


Figure  A.  7-3  Root  Locus  Compensator  Design 


A -35 


THE  ANALYTIC  SCIENCES  CORPORATION 


A.  3  SYSTEM  MODEL  SUMMARY 

In  Fig.  A.  8-1  the  complete  missile -target  intercept  model  is  portrayed 
with  all  of  the  subsystems  described  in  the  previous  sections  appropriately  inter¬ 
connected.  All  of  the  state  variables  are  depicted  except  Xg,  Xg,  and  (angle 
of  attack,  control  fin  deflection  and  the  seeker  compensation  state,  respectively) 
which  are  encompassed  in  the  linear  dynamics  represented  by  the  transfer  func¬ 
tions  g^(s),  g2(s),  and  (l+kg/s).  For  convenient  reference,  we  list  the  non- 
linearities  incorporated  in  the  system  model  and  indicate  their  form: 


Seeker  head  restoring  torque 


fl<*l)  =  k! 


®lim 


sign  (Xj) 


Seeker  gimbal  friction 


f2(x2)  =  k2  si^n  (x2> 


Receiver/signal  processing  characteristic 


V*>  ■ 


fUm  sign (t),  |' |  ><Um 


Range  dependent  noises  (2  nonlinearities) 

J1  2  /J2  2 

+  x.  _  Vx  +  y  +  x-  J  Vx  +  y 


Vs  ~  X16  ~r  "15 


14 
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iss4mbalancetorque -^3-nonIlnearities) 


Tm  =  mr0  trl  %  cos  («h  +  V  +  rl  (8m>2  si"  \  *  V 


+  vm  \  cos  -  9m  -  ®h  -  90>] 


LOS  angle  calculation 


Range  calculation 


0  =  tan”1  (y/x) 


J  2  2 

r  =  Vx  +  y 


Velocity  resolution  (4  nonlinearities) 

X  =  -  V  cos  (9  )  -  v,  cos  (9  ) 
m  it  a 

y  =  -  v  sin  (0. )  +  v.  sin  (0  ) 
J  m  4  t  a 


«  Acceleration  command  limiting 


ac  = 


a  sa 
1  c  max 


a  sign  (a  ')>  la '  I  > 

'  max  °  c  ’  c  max 


Proportional  guidance  law  with  secant  compensation 


c  =  "S  [vm  +vt 


COS  (0a  +  0) 


In  addition,  the  state  variables  are  defined  as  indicated  in  Table  A.  8-1. 
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TABLE  A.  8-1 

SYSTEM  STATE  VARIABLES 


State 

Variable 

Physical 

Notation 

Description 

xi 

9h 

Seeker  head  angle  with  respect  to  missile 
centerline 

X2 

• 

9h 

Seeker  head  angular  rate 

X3 

6 

m 

Missile  body  angle  with  respect  to  original 
LOS 

X4 

• 

8 

m 

Missile  body  angular  rate 

X5 

a 

Angle  of  attack  (incremental) 

X6 

6 

Missile  control  fin  deflection 

*7 

A 

8 

Estimated  LOS  angular  rate  (output  of  gui¬ 
dance  law  filter 

X8 

\ 

Missile  lead  angle  (angle  of  missile  velocity 
vector  with  respect  to  original  LOS) 

X9 

9a 

Target  aspect  angle  (angle  of  target  velocity 
vector  with  respect  to  original  LOS) 

x10 

y 

Cross -range  separation  (component  of  range 
normal  to  original  LOS) 

X11 

X 

Down-range  separation  (component  of  range 
along  original  LOS) 

x12 

ev 

Closing  velocity  error  state 

^3 

at 

Target  lateral  acceleration 

x14  tox16 

None 

Seeker  noise  states 

x17 

None 

Seeker  compensation  state 
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APPENDIX  B 

RANDOM  INPUT  DESCRIBING  FUNCTIONS  FOR 
THE  MISSILE -TARGET  INTERCEPT  PROBLEM 

B.  1  RANDOM  INPUT  DESCRIBING  FUNCTIONS  USED  IN  THE 
CADET  ANALYSIS 

The  random  input  describing  functions  (ridf's)  derived  in  this  sec¬ 
tion  are  based  on  the  assumption  that  all  system  variables  are  jointly  normal. 
To  summarize  the  quasi -linearization  procedure,  treated  more  fully  in 
Section  2.2  and  Appendix  C,  we  consider  a  general  function  of  two  variables*, 
cp(yi,  v2),  which  is  to  be  represented  approximate iy  by 

(p  (yi,  v2)  -  E[<p]  +  nr  (vj-m^+nj.  (v2-m2)  (B.l-1) 


Thus  is  it  necessary  to  evaluate 


The  nonlinearities  considered  are  generally  ordered  according  to  Appendix  A. 


In  this  appendix,  variables  are  generally  designated  v,  or  V]i  ,v2  etc.  as  re¬ 
quired,  to  avoid  the  need  to  refer  to  the  system  model.  The  statistics  are 
m  =  E[v],  p  =  E[(v-m)2],cr  =  vf!p  for  one  variable;  m^EfvJ,  pij  =  E[(vi-mi)(Vj-mj)], 
ai =  vPii  f°r  two  or  more  variables.  With  two  variables,  the  correlation 
coefficient  is  useful;  p  =  P12A/P11P22- 
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Missile-Target  Kinematics  —  The  nonlinearities  used  to  project 
the  missile  and  target  velocity  vectors  into  down-  and  cross-range  com¬ 
ponents  are  of  the  form  sin  (v)  and  cos  (v),  where  v  is  a  gaussian  random 
variable.  The  ridf's  required  for  these  functions  have  been  derived  pre¬ 
viously  to  be  (Refs.  3, 8) 


sin  (v) 


1  •  /  \  2 
=  —  sin  (m)  e 
m  v 

1 

"  2  P 

=  cos  (m)  e 


cos  (v)  - 


=  —  cos  (m)  e 
m 


=  -  sin  (m)  e 


(B.l-3) 


Guidance  Law  —  Referring  to  Eq.  (A.  5-4),  the  unconstrained  ac¬ 
celeration  command  is  a  highly  nonlinear  function  of  six  system  variables 
of  the  form 


v  •  Vi  +  V2  +  a3Vl 


(B.l-4) 


One  component  of  v'  is  a  product  of  states, 


=  v1v6 


(B.  1-5) 


and  a  second  constituent  is  given  by 


ccs  (Vo  +  8) 

<Mvi)  =  vi  - 7 - S7 

2  1  COS  (V3-  0) 


(B.  1-6) 
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when  8  is  given  by 


8  =  tan 


-1  /  V 


(B.  1-7) 


Substituting  Eq.  (B.l-7)  into  Eq.  (B.l-6)  yields 


VgCOS^)-  V4  sin  (V2) 
<P2(vi)  =  vi  v^cos  (v3)+v4sin(vq) 


(B.  1-8) 


For  the  product  of  variables,  <p^,  we  have  the  followinfc  results  (Ref.  3): 


E[Vv1’v6>]  =  nmiml+n: 


m„ 

m6  6 


=  mlm6  +  p16 


(1) 

\  "  me 

n^  =  m, 


(B.l-9) 


We  point  out  that  in  treating  nonlinear  functions  of  more  than  one  variable, 
it  is  generally  impossible  to  identify  the  individual  elements  of  the  describ¬ 
ing  function  for  the  mean  part,  e.  g. ,  in  this  case  we  do  not  obtain  nm^  and 
nm2  explicitly;  however,  the  quantity  E[oJ  is  all  that  is  required  for  propa¬ 
gation  of  the  mean  via  CADET  (refer  to  Eq.  (C.2-5)).  It  should  also  be  ob¬ 
served  that  the  representation  in  Eq.  (B.l-7)  is  not  useful  if  both  of  the 
means,  m^  and  mg,  and  the  covariance  are  zero  —  i.e.,  the  describing 
function  approximation  to  the  nonlinearity  is  zero  for  this  case.  Since  this 
point  is  considered  in  detail  in  Ref.  3,  we  do  not  treat  the  subject  here. 
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The  second  term^fri)  *n  Eq.  (B.l-4),  is  impossible  to  quasi- 
linearize  exactly  in  closed  form  under  the  gaussian  assumption;  thus  we  use 
a  generalization  of  the  series  approach  discussed  in  Section  5.2  (Ref.  9). 


E 


M = 


m. 

A-i  m.  1 
i  =  l  i 


5  5 


*2 

a  (p, 


,(2) 


n'  '  = 
r. 

l 


1  ys  y  s  ^2 

^(m.-m., . .  .,m_)  +  K- /J  /.  x — s — P-- 
2 '  1  2’  ’  5'  2  p,  am.  din.  i] 

d<pg ^ >^2*  *  *  *  ’  ^5^ 


(B.  1-10) 


9m. 


i  =  1,2, . . . , 5  (B. 1-10) 


Listing  the  partial  derivatives  called  for  in  Eq.  (B.  1-10)  requires  the  in¬ 
troduction  of  some  auxiliary  notation: 

^  =  m5  cos  mg  -  m^  sin  m2 


cos  m3  +  ra^  sin  mg 


il>3  =  -  mg  sin  nig  -  *n^  cos  m,. 


(B.  1-11) 


=  "m^  sin  m3  +  cos  m^ 


In  terms  of  these  expressions,  the  quantities  required  to  evaluate  Eq.  (B.l-10) 
can  be  shown  to  be 


<P2 (m1 , m2 , . .  .,m5) 


.  *l  ,  n(2) 

Sml  *2  rl 


(B. 1-12) 
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d<Po 


* 


3m, 


■y^m,  ^  n^ 

*2  1  r2 


3<0, 


*1*4 


3m, 


m<  -  n 


(2) 


Pc 


S«P0 


mimK 

- - sin  (m2  +  m3^  - 


3m, 


n 


(2) 


3<Po 


3mr 


mlm4 

2  sin  (m2+m3) 


—  n 


(2) 


-*2 

3  <p, 


3  m 


=  0 


»S  *3 


3m^3m2 


P 


2 

3  <Pf 


*1*4 


3m^3m3 


l/>r 


«.2 

3  <Pn 


mr 


3m^3m^ 


-  -~2  sin  (m2  +  m3) 
^2 


*2 

5  <Pn 


m. 


ra:  =  -^-sln(m2-m3) 


(B.  l-12)(cont. ) 
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amgdnig 


¥4 

~rm  1 
*2 


.2 

a  fj 
dm2Sm4 


mlm5 


COS  (m2  +mg) 


.2 

a  <P2  m1m4 

___  — 


2  5  £ 


cos  (mg  +m3) 


(B.  l-12)(cont. ) 


dm3am4  « 


mlm'i  /  \ 

=  . — (04  sin  (m2  +  m3)  - 


m,  m. 
1  4 


amoamc 
o  5 


$4  sin  (m2+  m3)  - 


6  *2  „  mlm5 


=  2  — g —  sin  (m3)  sin  (mg  +  m3) 


am.  am-  ,3 

4  5  02 


-  2m,.  cos  (m3)^  sin  (m2  +m3) 
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a2<p-  2m^m, 

- y  = - ^  cos  (m3)  sin  (m2+m3)  (B.l-12)(cont. ) 

3m5  *2 

Returning  to  Eq.  (B.l-4),  we  have 

E[<p7]  =  ^  E[pJ  +  agEf^j  +  a3m1  (B.l-13) 

where  E[c?-^]  and  Ej^]  are  given  in  Eqs.  (B.l-9)  and  (B.l-10);  the  random 
component  of  <p  can  be  expressed  in  terms  of  the  quasi-linear  gains  in  the 
same  equations  to  be 


(1)  (2) 

alnrx  +a2nr1+a3 

(2) 

a„n'  ’ 

2  r2 


a2nr 


a2n 


(2) 

3 

(2) 

*4 

(2) 


a2nr 


a<n^+  0  +  0 

1  r6 


=  rT  b  (B.  1-14) 

-  -<P 

Since  r^  is  a  quasi-linear  combination  of  the  random  components  of  the  six 
variables  o,  the  variance  is  approximately 
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p '  =  E  [r2  1  s:  E  rbT  rrTb  ]  = 
y<p  [_  J  |_-<P —  -<PJ 


P11  p12 
p12  P22 


p16  P26  ' 


.  .  p. 


(B. 1-15) 

Thus  we  can  obtain  the  approximate  mean  and  variance  of  the  unconstrained 
acceleration  command  in  terms  of  the  quasi-linear  gains  summarized  above. 
With  these  statistics,  we  treat  the  acceleration  command  limiter  in  the 
regular  manner;  i.  e. ,  we  assume  its  input  tp'  is  a  gaussian  random  variable 
with  statistics  m^  and  p^.  The  ridfs  for  the  limiter  are  given  by  (Refs.  2 


and  8) 

E[»]  -  n  m  -  f7  fcfeSV  G(/%^ 

m*  *  -/pn  v  yz  ;  \  & 


(B.  1-16) 


(<p..  +m'\  h..  -m'\ 

v  r^H-V)-1 

where  the  auxiliary  functions  G  and  PI  are  defined  in  terms  of  the  normal 
density  function,  PF,  as  follows: 


2 

Pm  =  -==  exp  (-4) 

you  '  1  / 


pm  =  f  pf(c) 


(B. 1-17) 


iji  PI(0)  +  PF(0) 
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The  describing  functions  outlined  above  complete  the  quasi- linear 
representation  of  the  guidance  law;  in  summary,  E[<p]  is  used  directly  in 
propagating  the  mean  vector  of  the  system,  and  the  random  component  of 
the  actual  acceleration  command  is 


r 

<P 


n 


r '  =  n 


r  o 
(p 


bTr 


(B. 1-18) 


and  nr^  b  ^  is  the  overall  describing  function  matrix  (in  this  case,  a  row 
vector)  specified  in  Eqs.  (B.  1-9),  (B.  1-10)  and  (B.  1-14). 

The  approach  outlined  above  in  Eqs.  (B.  1-4)  to  (B.1-*")  considers 
a  nonlinearity  of  the  form 

<p  =  <p(q'(vvv2,  •  ••  ,v6j) 


i.e. ,  a  nonlinear  function  of  a  nonlinearity.  Because  it  is  essentially  im¬ 
possible  to  quasi -linearize  this  relation  as  a  whole,  we  have  first  quasi- 
linearized  <p'  to  obtain  the  statistics  m'  and  p'  necessary  to  calculate  the 

ip  (D 

ridf  for  <p  (<p')>  Eq.  (B.l-16),  then  "cascaded  the  ridfs"  for  the  random  part 
in  arriving  at  Eq.  (B.  1-10).  While  this  is  not  a  completely  rigorous  pro¬ 
cedure,  we  rely  on  our  a  priori  knowledge  that  in  the  guidance  law,  <o'  is 
approximately  v^  times  the  closing  velocity  which,  while  not  a  constant,  is 
nearly  so  (i.e. ,  the  closing  velocity  has  a  mean  that  is  much  larger  than  its 
random  part  until  the  last  few  hundredths  of  a  second  of  an  engagement  in 
realistic  situations),  and  v^  is  generally  nearly  gaussian.  Thus  it  is  rea¬ 
sonable  to  assume  that  <p'  is  nearly  gaussian.  In  this  situation,  the  above 
technique  adequately  represents  the  guidance  law  nonlinear  effects,  as 
shown  in  Chapter  4. 

Range  Dependent  Noise  Sources  —  Since  the  range,  r  =  Jyp  +y2, 
enters  into  the  seeker  noise  model,  we  find  that  the  ridf’s  cannot  be 
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calculated  under  the  gaussian  assumption,  due  to  the  inability  to  evaluate 
certain  integrals  in  closed  form  (cf .  Section  5. 2).  The  first  nonlinearity 
we  treat  is  the  range  proportional  seeker  noise  term, 


(B. 1-19) 


where  we  assume  that  x  has  no  random  component,  as  before  (Section  5.2). 
The  expected  value  of  <o„,  given  by 


(B. 1-20) 


can  be  integrated  with  respect  to  v^  to  yield 


with  the  substitution  u  =  C^-mg)/^-  To  determine  the  ridf's  for  the  ran¬ 
dom  components,  we  take  the  indicated  partial  derivatives  with  respect  to 
m.;  these  are  simplified  significantly  by  the  further  assumption  that  the 
cross-range  separation,  y  =  V2,  has  a  zero  mean: 
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where  we  have  eliminated  terms  of  the  form 

_ 1  2 

ft  2  “2U  .  n 

\  uvto^u)  +m  e  du  =  0 
-go 


We  note  that  none  of  the  above  integrals  can  be  evaluated  in  closed  form. 
The  study  of  approximate  methods  presented  in  Section  5.2  demonstrates 
that  accurate  results  are  obtained  by  replacing  the  gaussian  pdf’s  in 
Eq.  (B.  1-22)  with  a  triangular  density  with  the  same  mean  and  standard 
deviation, 


|u|*/6 

|u|>/6 


(B. 1-21 


Furthermore,  we  recognize  that  nr^  is  the  expected  value  of  the  range,  as 
treated  in  Section  5.2,  so  we  immediately  have  (from  Eq.  (5.2-12)) 


m, 


Sk]--  sign  (m  ) 

°  V  A 


f~2  : 

V  l  +V  +  V 


+  ±LK{l+v2^ 

Je\  / 


(B.  1-24) 


nr  -  E[<p3]/m1 
where 


76  ct2 


(B. 1-25) 


Evaluating  the  last  integral  in  Eq.  (B.  1-22)  in  the  same  manner,  we  obtain 


A 


3  P12 


(1  +4y 


(B. 1-26) 
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The  inverse  range  proportional  noise, 


x  J 


J 


m?  + v? 
x  2 


can  be  treated  in  the  same  manner  as  <p_  above,  yielding 


fi+Ji7v2 


log (  |+l;-7l+? 


nr  —  E  [<p4]  /m3 


(B.l-27) 


n 


RecJyer/Signal  Processing  Characteristic  —  The  boresight  error 
limiter  has  the  same  form  as  the  acceleration  command  limiter;  conse¬ 
quently  the  evaluation  of  quasi-linear  gains  is  accomplished  with  Eq.  (B.l-16) 
substituting  the  boresight  error  statistics,  me  and  pc  for  m/  and  p^.  The 
boresight  error  is  a  linear  combination  of  the  LOS  angle,  missile  body  angle 
and  seeker  angle,  and  the  LOS  angle  is  related  to  the  system  variables  x  and  y 
by  the  inverse  tangent  relation,  Eq.  (A.  2-1).  Thus  to  determine  m€  and  pe, 
consider 


,  ,  -l/vl\ 

=  (V  '  V3  '  V4 

=  (al,v2)  -  V3-V4  (B.  1-28) 
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The  nonlinearity  is  quasi-linearized  below  so  that  the  mean  and  random 
components  of  the  boresight  error  are 


mf  =  E[<p5]  -  m3  -  m4 


r  = 


n<5>r1+n<Vr,-r4 
r^  1  r2  2  3  4 


(B.l-29) 


A  ■. 

"  —5  — 


The  boresight  error  variance  is  thus  approximately 


T 

P  -  bt 

-  5 


P11  P12 


14 


PI2  P22  •  '  p24 


P14  P24 


■44 


-5 


(B. 1-30) 


Using  mf  and  pf  from  Eqs.  (B.l-27)  and  (B.i-28),  the  boresight  error 
limiter  ridf's  can  be  evaluated  directly,  from  Eq.  (B.l-16). 


Inverse  Tangent  LOS  Angle  Calculation  —  The  nonlinearity 
<pp.  =  tan  *  cannot  be  quasi-linearized  in  closed  form  under  the 

gaussian  assumption;  thus  we  again  make  use  of  the  series  expansion  ap¬ 
proach  given  in  Lq.  i.B.1-10).  The  required  partial  derivatives  are 


^5 


2  2 
ml  +m2 


-m 


1 


2  ^  2 
m^  +  m2 


(B.l-31) 
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-2m^m2 

(mi +  mi)S 


ml-m2 


ml+m2 


(B.l-31)  (cont.) 


2m^m2 

R+ml) 


The  expected  value  of  0^  is  obtained  approximately  by  substituting  Eq.  (B.  1-31) 
into  the  series  expansion, 


E  [<p,-]  -  tan 


Isr)  +  ?V  -S\2  [mim2  (P22-Pll>+(ml-m2)I>12 

'  2/  (ml +  ”2)  J 


(B.  1-32) 


The  quantities  E[<pg],  nj^  and  n^  constitute  the  quasi-linear  approximation 
to  the  inverse  tangent  function. 


Seeker  Head  Restoring  Torque  —  The  nonlinear  representation  of 
the  "hard  spring”  effect  Qtudied  in  Section  4.4  is  of  the  form 


«>6  =  |v|ksign(v) 


(B.  1-33) 


where  k  is  an  integer.  The  ridf  is  directly  available  (Ref.  8),  in  general 
form  for  k  odd,  and  also  for  k  =  2,  as  follows: 


k  =  1 


n  =  1 

r 

n  =  1 

m 


(B.l-34) 
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i l 

3 

I 

H 


B 


!! 

S 

11 

n 

i 

4* 


sr* 
t  ' 

i-w 


A  *r 


k  =  2 


k  =  3 


k  =  5 


nr  =  4cr  PF  (“)  +  2m  ^2 PI  -  Ij 

%  ■  2oP1?(")+“>(1+©2)(2H(f)-1) 

n  =  3a2  +  3m2 

r  (B.l-34)  (cont.) 

—  Qrr^  4. 

n  =  ocr  +  m 
m 

4  2  2  4 

n  =  15a  +  30a  m  +  5m 
r 

4  2  2  4 

n  =  1 5a  +  10a  m  +  m 
m 


The  general  form  for  k  odd,  say  k  =  2g-f-l,  is 


n 


V  (2g+l)l  (1)(3)(5)  •  •  •  (2g-2j+l)  g-j  2j 
A  (2j)!  (2g  -  2j  +  1)!  p 


n 


m 


2g  v-  (2g+l)i  (1)(3)(5)-  •  •  (2g-  2j  -1)  g-j  2j 
m  *  fa  (2j+l)!  (2g-2j)!  P 


(B- 1-35) 


we  note  above  that  for  k  odd,  the  nonlinearity  in  Eq.  (B.l-33)  is  simply 
<pg  =  v2g+1  which  leads  to  quite  straightforward  polynomial  ridf's,  while  for 
k  even  the  factor  sign  (v)  leads  to  terms  involving  the  probability  function, 
PF,  and  probability  integral,  PI,  given  in  Eq.  (B.l-17). 


Seeker  Gimbal  Dry  Friction  —  The  dry  or  coulomb  friction  effect  is 
modeled  by  a  nonlinearity  of  the  form 


<p_(v)  =  sign  (v) 


(B. 1-36) 
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■which  is  often  used  to  represent  the  ideal  relay.  The  ridf's  for  this  effect 
are  available  in  Ref.  8, 


n 


n 

m 


1_ 

m 


(B.l-37) 


Seeker  Mass  Imbalance  —  The  effective  disturbance  torque  due  to 
an  offset  of  the  seeker  center  of  gravity  from  the  gimbal  has  three  terms, 
with  two  distinct  forms: 


(P8  =  cos  (v2) 

2 

<P9  =  vx  sin  (v2)  (B.  1-38) 

In  the  nonlinearities  <pg  and  <pg,  the  variables  v^  and  v2  are  generally  not 
state  variables,  but  linear  combinations  thereof,  viz., 


v  = 

- 1 

< 

_ _ 1 

r 

XI 

1 

1 _ 

+ 

0 

v2 

- 

-2  — 

l 

a> 

o 

i - 

(B. 1-39) 


where  we  have  indicated  that  the  angle  v2  generally  contains  a  constant  param¬ 
eter,  0q,  which  is  the  angular  coordinate  of  the  seeker  center  of  gravity  with 
respect  to  the  antenna  centerline  (Fig.  A. 6-1).  Then  the  statistics  of  v.^  and 
v2  can  be  obtained  from  the  state  variable  statistics,  m  =  E[x]  and 
P  =  E  [(x-m  )(x-m  )T],  as 


m 


v 


0 


H  m  ± 


(B.l-40) 
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where 


and 


(B.l-41) 


P 

v 


P11  P12 
P12  p22 


T 

HPH 


(B. 1-42) 


The  first  nonlinearity  (and  its  companion  form,  <p =  v^sin  v2) 
have  been  treated  in  Ref.  3  by  considering  the  real  and  imaginary  part  of 
the  function  exp  (jv2),  as  follows:  define 


*  A 

*8  =  V 


=  Vj  (cos  v2  +  j  Sin  v2)  =  +  j(f>. 


10 


(B.  1-43) 


The  expected  value  of  this  complex  function  can  be  evaluated  to  be  (Ref.  11) 


E 


;«$]- 


P22 

2 


(mx+  jp]2)  (cos  m2+j  sin  m2)  (B.l-44) 


From  this  expression,  the  required  ridis  are  available  by  inspection: 


E  [<pg  ]  =  Re 


(■KD- 


22 


(raj  cos  m2  -  pj2  sin  m2)  (B.l-45) 
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'22 


n  =  e  cos 


m2 


■22 


(B.l-45)(cont.) 


n 


-  e  (m^  sin  mg  +  p^2  cos  m2) 


/  \  22 

E[<pi0]  *  ImfE^  ])  =  e  2  (“l  sin  m2  +p12  cos  m2) 


* 22 

„  _  "  2  . 

n  -  e  sin  m0 

r  1  1 


(B.l-46) 


22 


n 


=  e  (m^  cos  m2  -  p12  sin  m2) 


The  same  approach  is  utilized  in  evaluating  ridf's  for  <pQ;  defining  its  com- 

2  y 
panion  form  to  be  <p^  =  cos  Vg  yields 

*  2  ^v2  2 

ip^  =  vi  e  =  (cos  v2  +  j  sin  v2) 


yields 


=  <Pn-M>9 


p22 

2 


(B.  1-47) 


=  e  2  (mi  +P1l“P12)  +  2imlp12  (cos  m2 +  J sin  m2)  (B.l-48) 


/ 
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Thus  the  quasi-linear  representation  of  <pg  and  <p^  require 


E  [<p11]  =  Re  (e 


8 


=  e 


■22 

2 


[(ml  +  P).l-p12)C0S  m2  -2mlp12  sln  “2] 


•22 


n 


=  e  (2m^  cos  m2  -  2pj^  s*n  m2)  (B.l-49) 


22 


n  =  -e 


x2 


(m^  +  pn  -  p12)  sin  m2  +2m1p12  cos  m2 j 


P 


22 

2 


(^  \  -  9  '  t  o  \ 

E  <pg  j  =  e  2miPi2  cos  m2  +  \mi  +pll  "P12/Sin  m2 


■22 


n  =  e 


|^2p^2  cos  m2  +2m^  sin  m2J  (B.l-50) 


n  =  e 


2 


_22 

2  -2m1P12  sin  m2  +(mi  +  pxx  -p12)  cos 


This  completes  the  exposition  of  the  ridf's  necessary  for  the  statistical 
analysis  of  the  missile -target  Intercept  model  treated  in  the  present  study. 
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B.2  RANDOM  INPUT  DESCRIBING  FUNCTIONS  USED  IN 
SENSITIVITY  ANALYSIS 


In  Section  5.2  we  consider  three  nonlinearities,  the  limiter,  power 
law  nonlinearity,  and  sinusoidal  operator,  and  portray  the  variation  caused 
in  ridf  calculation  by  varying  the  probability  density  function  (pdf)  of  the 
input.  Seven  quite  different  pdfs  were  studied,  denoted  p.(x),  i=  1, 2, . . . ,  7, 
as  given  in  Eq.  (5.1-4);  each  pdf  leads  to  a  describing  function  given  by 


(B.2-1) 


where  x  is  a  zero-mean  random  variable  with  standard  deviation  a. 


Limiter  —  For  this  nonlinearity,  all  ridfs  are  a  function  of  the 
ratio  of  the  input  standard  deviation  to  the  break  point  (Fig.  5.1-2),  here 
denoted  by 

H  =  c/6 

In  terms  of  this  parameter,  we  have 


n 


r 


3 


1, 

1  (*/6  m  +  1) 

1 "  4 

36  ji 


C/6  p-1), 


-}=<»<«> 


(exponential) 

(gaussian) 


(triangular) 
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Random  input  describing  functions  are  quite  complicated  when  both  the  pdf 
and  nonlinearity  are  comprised  of  several  linear  segments,  as  is  the  case 
for  p4,  p6  and  p^. 

Power  Law  Nonlinearity  —  The  nonlinearity 

2 

f(x)  =  x  sign  (x) 

leads  to  particularly  convenient  ridf’s  for  the  pdf’s  chosen  in  this  study. 


m)Im 

Q 

II 

2.12  a 

(exponential) 

4 

T  or  = 
-/2tt 

=  i.60cr 

(gaussian) 

376  _ 

—  a  ’ 

1.47  a 

(triangular) 

= 

575 

1.39  a 

(P4) 

. 

4  a  - 

1.30  a 

(uniform) 

= 

777 

1.19  a 

(Pfc) 

54  76  „  . 
125 

=  1.06  a 

(P7) 

Sinusoidal  Operator  -  For  the  nonlinearity 

f(x)  =  sin  x 

the  following  ridf’s  are  obtained: 
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n  = 

(exponential) 

rl 

2 

a 

3 

CO 

II 

e  2 

(gaussian) 

n  = 

r 

2(1  -cosx3)-x3  sin x3  ; 

x3  =  0 

(triangular) 

3 

3a  L  J 

II 

j^2(cosx4  "  cos  (3x4))+  x4(sinx'4  -  3sin(3x4))j  ; 

x.  =  Vo. 6  a 

4 

(P4) 

n  = 
r5 

if.  2  I 

— -r  x.  sin  xc  -  x_  cos  x-  : 

3a4  L  5  55  5J 

xc  =  Jia 

D 

(uniform) 

n  = 
r6 

-^4  x  (sinx6-sin(2x6))-(l-2cosx6  +  co£:(2x6))  ; 
3a  -  J 

(P6) 

n  = 
r7 

9^  r  0  0 

-=•4  x  (2sin(2x7)-  |Sinx7  -jsintfxrj)) 

3a  L 

X 

-O 

II 

H^) 

Q 

+2(cos  (2x7)-  cosx7  -  cos  (3x7))j  ; 

(P7) 

The  ridf'a  for  the  three  nonlinearities  considered  here  are  plotted  in 
Figs.  5.1-2  (limiter),  5.1-4  (power  law),  and  5.1-6  (sinusoidal  operator). 


B-23 


THE  ANALYTIC  SCIENCES  CORPORATION 


APPENDIX  C 

THE  COVARIANCE  ANALYSIS  DESCRIBING 
FUNCTION  TECHNIQUE  (CADET) 


The  Covariance  Analysis  Describing  function  Technique  (CADET)  is  a 
method  for  the  direct  determination  of  the  statistical  properties  of  a  nonlinear 
system  with  random  inputs,  recently  conceived  and  developed  at  The  Analytic 
Sciences  Corporation  (Refs.  1,  2,  and  3).  The  principal  advantage  of  this  tech¬ 
nique  is  that  it  obviates  monte  carlo  simulations,  thereby  achievin''  substantial 
savings  in  computer  processing  time.  We  first  motivate  the  discussion  by  re¬ 
viewing  the  covariance  analysis  method  for  linear  svstems;  then  we  develop  an 
analogous  procedure  (CADET)  for  the  nonlinear  case. 

C.  1  COVARIANCE  ANALYSIS  FOR  LINEAR  SYSTEMS 

The  dynamics  of  a  linear  continuous -time  stochastic  system  can  be  re¬ 
presented  by  a  first-order  vector  differential  equation  in  which  x(t)  is  the  system 
state  vector  and  w(t)  is  a  random  forcing  function  vector, 

x(t)  =  F(t)  x(t)  +  G(t)  w(t)  (C.  1-1) 


Figure  C.  1-1  illustrates  the  equation.  The  state  vector  is  composed  of  any  set 
of  variables  sufficient  to  completely  describe  the  behavior  of  the  system.  The 
forcing  inaction  vector  w(t)  represents  disturbances  as  well  as  control  inputs  that 
may  act  upon  the  system.  In  what  follows  the  forcing  function  w(t)  is  assumed  to 
be  composed  of  a  mean  value  b  and  random  component  u,  .he  latter  being  com¬ 
prised  of  elt.nentf  which  are  uncorrelated  in  time;  that  is,  u(t)  is  white”  noise 
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Figure  C.  1-1  Continuous  Representation  of  the  Linear 
Dynamic  System  Equations 

having  the  spectral  density  matrix  Q(t);  thus  w(t)  is  specified  by* 


E[w(t)]  =  b(t) 

E[u(t)  uT(r)]  =  Q(t)  6  (t-r) 


(C. 1-2) 


Similarly,  the  state  vector  has  a  mean  component  m(t),  and  a  random  part 
r  =  x  -  m  ;  x,  then,  may  be  described  statistically  by  its  mean  and  covariance 
matrix, 


m(t)  =  E[x(t)] 

P(t)  *  E[r(t)  rT(t)] 


(C. 1-3) 


The  equation  for  the  propagation  of  the  mean  vector  and  covariance 
matrix  for  the  system  described  by  Eq.  (C.  1-1)  can  be  written  as  (Ref.  7). 


* 

E  denotes  ensemble  expectation,  or  average  value;  T  denotes  the  matrix 
transpose. 
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m  =  F(t)  m  +  G(t)  b 


P(t)  =  F(t)  P(t)  +  P(t)  F(t)T  +  G(t)QGT(t) 


(C. 1-4) 


The  elements  of  m(t)  represent  the  effects  of  deterministic  initial  conditions 
and  biases  due  to  the  system  inputs  (b£0).  The  diagonal  elements  of  P(t)  are 
the  mean  square  values  of  the  random  components  of  the  state  variables;  the 
off-diagonal  elements  represent  the  degree  of  correlation  between  the  random 
components  of  the  various  state  variables. 

Equation  (C.  1-4)  provides  a  direct  method  for  analyzing  the  statistical 
properties  of  x(t).  This  is  to  be  contrasted  with  the  monte  carlo  method,  where 
many  sample  trajectories  of  x(t)  are  calculated  from  computer-generated  ran¬ 
dom  noise,  or  random  numbers  in  the  case  of  a  digital  computer,  using  Eq. 

(C.  1-1).  If  q  such  trajectories  are  generated  —  denoted  by  xk(t),  k  =  1, . . . ,  q 
--  then  m(t)  and  P(t)  are  given  approximately  by 

1  q 

m(t)  =  m(t)  =  -  2  xk(t) 

4  k=l 

(C.  1-5) 

P{t)  =  P(t)  =  q^y  £  r k(t)  rkT(t) 

k  =  l 

where  r  k  =  xk  -  m.  In  the  limit  as  q  approaches  infinity,  we  are  assured  that 


lim  m(t)  =  m(t) 
q>co 

lim  P(t)  =  P(t) 

q->co 


(C. 1-6) 
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Note  that  Eq.  (C.  1-4)  provides  exact  solutions  for  m(t)  and  P(t),  to 
within  computer  integration  accuracy,  whereas  the  monte  carlo  method  yields 
approximate  solutions  for  a  finite  value  of  q.  Furthermore,  Eq.  (C.  1-4)  need 
be  solved  only  once  over  the  time  interval  of  interest,  whereas  Eq.  (C.  1-1) 
must  be  solved  repeatedly  using  the  monte  carlo  technique;  consequently  the 
direct  analytical  method  is  not  only  exact,  but  is  also  generally  the  most  effi¬ 
cient  technique  for  analyzing  linear  systems.  Our  purpose  here  is  to  describe 
a  procedure  whereby  the  statistics  of  a  nonlinear  system  can  be  computed  ap¬ 
proximately  using  a  recursion  relationship  similar  in  form  to  Eq.  (C.  1-4). 

C.  2  COVARIANCE  ANALYSIS  FOR  NONLINEAR  SYSTEMS 
The  nonlinear  counterpart  of  Eq.  (C.  1-1)  is 

x(t)  =  f(x(t))  +  Gw(t)  (C.  2-1) 


as  depicted  in  Fig.  C.  2-1.  In  order  to  develop  a  covariance  analysis  method 
similar  to  that  used  for  linear  systems,  it  is  desirable  to  approximate  f(x)  in 
Eq.  (C.  2-1)  as  a  linear  operation  on  x(t).  The  success  of  CADET  in  achieving 
its  goals  --  providing  an  efficient  and  accurate  means  for  the  direct  statistical 
analysis  of  system  performance  —  depends  on  how  well  f(x)  can  be  approximated. 
Thus,  consider  approximating  the  nonlinear  function  J[(x)  in  Fig.  C.  2-1  by  a 
linear  function,  in  the  sense  suggested  by  Fig.  C.  2-2.  The  approximating  out¬ 
put,  f  (x),  is  comprised  of  the  sum  of  two  terms,  one  linearly  related  to  m  and 
the  other  linearly  related  to  r.  The  gain  matrices,  Nm  and  Nf,  are  chosen  to 
minimize  the  generalized  mean-square  error  between  the  actual  and  approxi- 
mate  outputs  defined  by  E[e  SeJ  when  S  is  any  symmetric  positive  definite  ma¬ 
trix.  This  is  often  referred  to  as  a  quasi-linear  approximation  technique  and 
Nm  and  N^  are  called  generalized  describing  function  matrices. 
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Calculation  of  N  and  N  is  readily  accomplished.  Note  first  from 
m  r  " 

Fig.  C.  2-2  that 

e  =  f(x)  -  f  Jx) 


=  f(x)  -  N  m  -  N  r 
—  —  m —  r- 


(C. 2-2) 


T 

Formally,  it  can  be  shown  that  E  [e  S  e]  is  minimized  if  Nm  and  Nr  satisfy 
(Ref.  1) 


Nmmm^  =  E [  _f(x) ] 
NE[rrT]  =  E[f(x)rT] 


(C.  2-3) 


The  second  relation  directly  defines  Nr, 

N.  =  E[f(x)rT]  P'1 


(C.  2-4) 


since  a  unique  P~*  always  exists.  Rather  than  attempting  to  solve  for  N  -- 

T  ^ 

which  requires  a  pseudo  inverse  since  (mm  )  is  always  singular  --  simply  note 
that 


NmH  =1  =  E[f(x)j 


(C.  2-5) 


This  result  is  all  that  is  required  to  solve  the  problem  at  hand,  as  shown  below. 

Evaluatingthe  expectations  in  Eqs.  (C.  2-4)  and  (C.  2-5)  requires  an  assump¬ 
tion  about  the  probability  density  function  of  x(t).  Most  often  a  gaussian  densi*',r 
is  assumed,  although  this  need  not  be  the  case  from  a  theoretical  point  r 
However,  the  gaussian  assumption  leads  to  the  significant  computatio  pli- 

fication  discussed  in  Section  C.  4.  In  addition,  it  is  often  physically  reasonable 
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to  assume  that  x  is  gaussian,  because  the  superposition  effect  created  by  the 
linear  part  of  the  system  dynamics  tends  to  transform  the  nongaussian  outputs 
of  the  nonlinearities  into  approximately  gaussian  state  variables,  by  the  central 
limit  theorem  (Ref.  9). 

C.  3  DEVELOPMENT  OF  THE  MEAN  AND  COVARIANCE  EQUATIONS 

Replacing  the  nonlinear  function  of  Fig.  C.  2-1  by  the  describing  func¬ 
tion  approximation  given  in  Fig.  C.  2-2  results  in  the  quasi-linear  system  model 
illustrated  in  Fig.  C.  3-1.  The  quasi-linear  system  model  consists  of  two  parts; 
one  part  to  propagate  the  mean  or  deterministic  portion  of  the  signal,  and  the 
other  to  propagate  the  random  portion.  Due  to  the  nonlinearity  of  the  system, 
the  two  loops  are  coupled  through  the  describing  function  matrices;  under  the 
gaussian  assumption,  both  matrices  are  functions  only  of  the  mean  and  covar - 
fence.  The  differential  equation  for  the  propagation  of  the  mean  vector  for  the 
quasi-linear  system  of  Fig.  C.  3-1  is 

m  =  Nm  P)  m  +  Gb  =  i  +  Gb  (C.  3-1) 


R -11951 


Figure  C.  3-1  Quasi-Linear  System  Model 
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The  system  covariance  propagation  equation  is  obtjiined  by  computing  the  covar¬ 
iance  of  the  random  signal  component,  r,  using  linear  system  covariance  analy¬ 
sis  as  discussed  in  Section  C.  1.  The  resulting  matrix  covariance  equation  is 

P  =  Nr(m,P)  P  +  PNrT(m,  P)  +  GQGT  (C.  3-2) 


These  equations  are  nonlinear  and  coupled  through  the  mean  and  random  com¬ 
ponent  describing  function  matrices,  N  and  N  ,  respectively.  Initialization 

m  r 

requires  only  the  initial  condition  statistics, 


m(0)  =  E[x(0)] 

P(0)  =  E[r(0)rT(0)] 


(C. 3-3) 


Equations  (C.  3-1)  and  (C.  3-2)  together  with  Eqs,  (C.  2-4)  and  (C.  2-5)  are  the 
key  equations  of  CADET. 


C.  4  SPECIAL  RELATIONSHIPS 

A  few  special  relationships  are  worth  noting.  First,  we  observe  that 
when  the  system  is  Unear  [f(x)  =  Fx],  Eqs.  (C.  2-4)  and  (C.  2-5)  immediately 
lead  to  the  result  Nm  =  N  =  F.  Hence,  Eqs.  (C.  3-1)  and  (C.  3-2)  reduce  to  the 
familiar  equations  for  the  propagation  of  the  mean  and  covariance  in  linear  sys¬ 
tems,  indicated  in  Eq.  (C.  1-4).  For  a  nonlinear  system,  computation  of  the 
describing  function  matrices  in  Eqs.  (C.  2-4)  and  (C.  2-5)  requires  an  assump¬ 
tion  about  the  form  of  the  probability  density  of  the  input  signal  to  the  nonlinear¬ 
ity.  This,  it  should  be  noted,  is  the  only  assumption  required  in  the  CADET 
formulation.  The  assumption  of  a  jointly  gaussian  probability  density  funetioi 
for  the  elements  of  x  is  most  attractive  from  two  viewpoints.  First,  the  cal¬ 
culation  of  Nm  and  Nr  tends  to  be  simplest  under  this  assumption.  In  addition, 
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many  dynamic  system  models  tend  to  contain  more  linear  than  nonlinear  ele¬ 
ments.  This  is  advantageous,  since  low-pass  linear  filtering  is  necessary  to 
insure  that  nongaussian  nonlinearity  outputs  result  in  nearly  gaussian  nonlinear¬ 
ity  inputs  as  signals  propagate  through  the  system.  This  so-called  "filter 
hypothesis"  is  common  to  all  describing  functions  analyses  (Ref.  8). 

One  result  of  the  gaussian  assumption  is  that  can  be  directly  com¬ 
puted  (Ref.  16)  from  the  relationship 

Nr(='  P)  =  4  E[I<-)]  =  4  [N«>-]  (C'  4"3) 

This  is  indeed  a  useful  relationship,  since  calculating  Nmm  is  required  in  any 
case  for  the  propagation  of  the  mean  (Eq.  (C.  3-1)).  It  is,  in  practice,  much 
easier  to  employ  Eq.  (C.  4-3)  than  to  solve  Eq.  (C.  2-4)  for  Nf(m,  P).  A  direct 
result  of  Eq.  (C.  4-3)  is  the  fact  that  Nm  and  Nr  can  be  formed  by  first  replacing 
the  individual  nonlinear  elements  of  f(x)  by  the  appropriate  scalar  describing 
function  gains.  This  is  an  extremely  powerful  property  since  a  large  number  of 
describing  functions  have  been  catalogued  in  Ref.  8;  consequently  and 
can  be  formed  directly  from  inspection  of  the  system  equations  or  block  diagram, 
in  many  cases. 
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